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ABSTRACT 


An  updated  version  of  an  earlier  waveguide  program  has  been  vrltten  In 
the  FORTRAN  compiler  language.  The  new  program  Includes  the  following 
alterations  or  additions;  (1)  An  Improved  Runge-Kutta  routine;  (2)  A 
procedure  which  finds  directly  the  starting  solutions  at  the  top  of  the 
Ionosphere;  (3)  Provision  for  Including  up  to  five  charged  particles;  • 

(4)  A  method  for  economically  finding  approximate  waveguide  solutions; 
and  modifications  for  use  of  the  program  for  ELF  wave  propagation. 
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I.  INTRODUCTION: 

This  report  is  the  latest  in  a  series  (1,  2,  3,  4,  5,  7)  which  de¬ 

scribe  computer  programs  which  have  been  developed  for  calculating  iono¬ 
spheric  reflection  coefficients  and/or  waveguide  modal  parameters.  All 
earlier  computer  programs  were  written  in  the  NELIAC  computer  compiler  lan¬ 
guage.  This  report  describes  a  program  written  in  FORTRAN  63  for  the  CDC 
computer  which  solves  the  modal  equation  for  propagation  within  an  earth 
ionosphere  waveguide.  The  experienced  programmer,  familiar  with  the  FORTRAN 
language,  should  find  it  easy  to  make  the  necessary  modifications  for  use  of 
the  program  in  other  computers. 

The  solutions  of  the  modal  equation  are  determined  from  the  reflection 
coefficient  matrixes  of  the  ionosphere  and  of  the  earth  as  seen  looking  up¬ 
ward  and  downward  from  a  position  within  the  waveguide.  The  program  system¬ 
atically  solves  the  modal  equation  for  as  many  modes  as  are  desired.  The 
propagation  parameters,  attenuation  rate,  phase  velocity,  and  excitation 
factor  are  determined  for  each  solution  of  the  modal  equation.  The  method 
makes  full  allowance  for  both  earth  curvature  in  the  direction  of  propagation 
and  the  orientation  and  intensity  of  the  geomagnetic  field  with  respect  to 
the  path  of  propagation. 

The  reflection  coefficient  matrix  of  the  ionosphere  is  obtained  numerically 
following  methods  developed  by  Budden  (8).  The  ionosphere  is  described  in 
terms  of  the  vertical  distribution  of  charged  particle  density,  the  approp¬ 
riate  charged  particle  neutral  particle  collision  frequencies  as  a  function  of 
height,  and  the  Intensity,  dip  angle,  and  azimuth  of  the  geomagnetic  field. 

As  many  as  five  charged  particles  of  either  sign  and  of  arbitrary  mass  may  be 
Included  as  constituents.  At  the  option  of  the  user  the  program  may  be  used 
to  calculate  and  print  a  listing  of  the  ionospheric  reflection  coefficients  only. 

In  addition  to  the  real  constituents  the  medium  contains  a  fictitious 
linear  distribution  of  permittivity  extending  to  the  ground  which  simulates 
(takes  into  account)  the  effects  of  earth  curvature. 

The  equations  used  for  obtaining  the  reflection  coefficient  matrix  looking 
at  the  ground  depend  on  the  height  at  which  the  modal  equation  is  solved.  If 
the  modal  equation  is  evaluated  below  the  ionosphere  but  above  the  earth,  be¬ 
cause  of  the  linear  distribution  of  permittivity,  the  surface  reflection 
coefficients  are  expressible  in  terms  of  solutions  to  Stokes’  equation  and 


their  derivatives.  If  the  modal  equation  is  evaluated  at  the  ground  all 
earth  curvature  effects  are  included  in  the  ionospheric  reflection  coeffici¬ 
ents  and  ordinary  Fresnel  reflection  coefficients  are  obtained  at  the  ground. 
The  ground  conductivity  is  an  arbitrary  input  so  that  results  represent  pro¬ 
pagation  over  a  realistic  finitely  conducting  earth. 

The  modal  equation  is  especially  useful  for  calculating  the  effect  on 
elf,  vlf ,  and  If  guided  propagation  by  anomalous  conditions  such  as  caused 
by  solar  flares  or  nuclear  explosions.  Step  models  or  exponential  models  of 
the  ionosphere  are  often  so  unrealistic  as  to  be  virtually  useless  for  eval¬ 
uating  the  effects  of  such  anomalies.  The  present  program  calculates  the 
attenuation,  phase  velocity,  and  excitation  of  Important  waveguide  modes  for 
realistic  waveguide  parameters.  The  program  thus  allows  one  to  make  realistic 
estimates  of  changes  in  field  strength  and  phase  velocity  over  critical  paths. 

In  susmary,  the  physical  Input  parameters  for  the  waveguide  program  are: 

1.  The  azimuth  of  propagation  with  respect  to  the  horizontal  component 
of  the  geomagnetic  field. 

2.  The  geomagnetic  dip  angle  and  field  strength. 

3«  The  radio  frequencies. 

4.  Vertical  profiles  of  up  to  five  charged  species  (electron,  positive 
and  negative  ions). 

5*  An  exponential  but  otherwise  arbitrary  collision  frequency  distribution 
with  height  for  each  specie  of  charged  particles. 

6.  The  ground  conductivity  and  dielectric  constant. 

The  output  can  be  either  a  printout  of  the  reflection  coefficient  matrix 
of  the  ionosphere  at  desired  height  increments  within  and  below  the  ionosphere 
or,  by  using  the  complete  program,  the  waveguide  mode  parameters  for  as  many 
modes  and  radio  frequencies  as  desired. 


II.  THE  IONOSPHERE  REFLECTION  MATRIX 
A.  Exact  Solution 

The  Ionosphere  reflection  matrix,  R  ,  at  height,  d.  Is  obtained  by  a 
numerical  Integration  of  differential  equations  given  by  Budden  '  .  The 
differential  equations  are  of  the  form 
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Because  of  some  symmetry  In  the  elements  of  the  S -matrix.  It  may  be  written 
in  the  form 
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where  C  =  cos  6  and  Q  Is  the  complex  angle  of  Incidence.  The  T -matrix  Is 
of  the  form 
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6e  c  electric  permittivity  of  free  space 
yClg  -  magnetic  permeability  of  free  space 
Jy  =  geomagnetic  field 

CO  =  angular  wave  frequency 
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The  quantities  N,  e,  m,  (nu)  have  different  values,  at  a  given  height, 
for  each  species.  The  quantities  are 


N  ■  species  density 

e  =  charge  of  electron  or  of  ion  (negative  for  electrons  and  negative  ions) 
m  =  mass  of  electron  or  of  ion 
collision  frequency 

For  each  species  there  will  be  an  jM^.  The  matrix,  _M,  used  in  computing 
the  T -matrix  is  then 


(II-5) 


The  differential  equations  are  integrated  by  the  Runge-Kutta  technique, 
starting  at  some  height  above  vhich  negligeble  reflection  is  assumed  to  take 
place.  Error  control  is  by  means  of  a  comparison,  for  each  step,  of  the 
increments  in  the  elements  of  ^  computed  with  the  fourth  order  Runge-Kutta 
method,  and  that  computed  vlth  a  second-order  integration  step.  Although  not 
a  completely  rigorous  error  test,  it  has  been  found  to  be  very  reliable. 

The  eight  integration  variables  are  the  arguments  and  the  logs  of  the 
magnitudes  of  each  of  the  four  elements  of  the  reflection  matrix.  The 
derivatives  of  these  variables  cure  obtained  vlth  the  expression 


(Il-t>) 


implies  the  complex  logarithm  of  the 

It  turns  out  that  the  initial  conditions  for  the  integration,  l.e.,  the 

starting  value  of  &  is  the  value  of  for  a  sharply-bounded  ionosphere  vlth 

properties  corresponding  to  the  ionosphere  at  the  top  of  the  given  electron 

density  and  collision  frequency  profiles.  The  solution  used  is  that  described 

(9) 

in  a  paper  by  C.H.  Sheddy,  '  '. 
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B.  Approximate  Solution 

Most  of  the  running  time  of  the  NELC  waveguide  program  is  in  the  compu¬ 
tation  of  £  which  is  used  in  the  modal  equation  |£&-i|=o 
Since  the  eigen  values  of  9  usually  occur  along  a  line  in  the  complex  e-plane, 
it  often  suffices  to  compute  ^R,  for  two  or  three  values  of  9  and  interpolate 
between  these  values  to  obtain^?  -  f(e).  In  this  way  the  total  nan  time  is 
cut  to  a  fraction  of  that  for  the  exact  solution. 

Approximate  solutions  may  be  used  for  exploratory  work  where  one  does 
not  require  great  accuracy.  It  may  also  be  used  in  cases  where  the  electron 
density  profile  itself  is  not  known  too  well,  since  the  errors  resulting  from 
use  of  the  interpolation  may  be  no  greater  than  the  uncertainties  due  to  not 
knowing  the  profile  too  well.  Even  when  exact  solutions  are  desired,  approxi¬ 
mate  solutions  may  be  obtained  first  as  starting  values  for  exact  solutions. 

The  Lagrange  interpolation  polynomial  is  used  in  the  program: 
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where  log  R  implies  the  complex  log  of  an  element  of  ^R.  is  a  value  computed 
with  the  full-wave  solution  at  the  angle,  9^. 

The  user  must  select  n  values  of  9^  when  using  the  approximate  scheme. 

Three  or  four  values  of  9 k  generally  give  best  results,  often  two  suffice. 

For  best  results  these  should  be  chosen  as  being  in  a  line  in  the  complex  9- 
plane  roughly  approximating  the  line  along  which  the  eigen  solutions  are  ex¬ 
pected  to  be.  It  also  has  seemed  best  not  to  have  more  than  one  9k  with  the 
same  value  for  the  real  parts  or  the  same  value  for  the  imaginary  parts. 

When  using  the  waveguide  program  for  the  first  time,  the  user  might  try 
several  sets  of  0^  for  the  same  set  of  input  data  in  order  to  gain  a  feel  for 
the  sensitivity  of  the  approximate  eigen  solutions  on  the  9^. 
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III.  SUMMARY  OF  MODE  PARAMETERS 

Following  the  formalism  developed  by  Sudden,  whereby  one  wwks  within 
the  framework  of  a  planar  stratified  medium  and  includes  earth  curvature  via 
the  artifice  of  a  modified  refractive  index,  it  is  well  known  that  the  mode 
equation  may  be  written  as  the  deteroinantal  equation. 


FD(e)  =  I R^©)  Rp  (©)  -  l  =0 
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where 
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is  the  reflection  matrix  looking  up  into  the  ionosphere  from  height  D  and 
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is  the  reflection  matrix  looking  down  from  height  D  towards  the  ground.  The 
angle,  0,  is  the  angle  of  incidence  at  height  z  =  H  where  the  modified  index 
of  refraction 
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is  unity.  The  full  wave  equations  programed  for  the  purpose  of  determining  the 
elements  of  Rp(©)  are  given  in  section  II.  It  is  emphasized  that  earth  curva¬ 
ture  is  included  in  the  calculations  of  the  elements  V.M ,  by  the  suscepti¬ 
bility  matrix.  Since  is  assumed  to  be  diagonal  it  is  essential  that  the 

level  D  be  chosen  at  or  below  the  base  of  the  ionosphere.  Except  for  D  -  0  and 


x 


some  ELF  calculations,  Rjj,  is  calculated  as  described  in  reference  ^  in 

terms  of  solutions  to  Stokes'  equation  and  their  derivatives.  The  solution  to 

Stokes'  equations  which  we  have  chosen  to  work  with  are  the  modified  Hankel 

12 

functions  of  order  1/3  as  defined  in  reference  .  It  should  be  observed  that 
the  latter  functions  are  linearly  related  to  the  perhaps  better  known  Airy 
functions.  When  D  =  0  or  when  the  magnitude  oi  the  hyperbolic  sine  of  the 
imaginary  part  of  the  eigen  angle  exceeds  a  preassigned  limit,  the  matri> 

is  calculated  directly  from  the  Fresnel  reflection  formulas. 

The  principal  objective  of  the  present  work  is  the  solution  of  (1)  for  the 
eigen  angles,  dn«  To  achieve  this  the  well  known  Newton's  method  is  used.  The 
procedure  is  to  guess  a  solution  0.  to  the  equation  f^»(©)  =  0.  The  function  Fd(& 
is  then  reevaluated  for  0.  +  §0  and  the  correction  to  0,  found  from  the  equation 

A©  =  -  _ F(Q)  ^0 _  (III-5) 

F(©o+  &0)  -  F(©0) 

The  correction  determined  by  Eq.  (5)  is  then  evaluated  and  the  process  re¬ 
peated  until  the  quantities  (A0r|  and  |  A©i|  are  reduced  to  within  a 
preassigned  tolerance.  For  the  first  iteration  the  value  of  ^0  is  taken 
equal  to  0.01*  •  The  subsequent  choices  of  ^0  are  contigent  upon  proximity  of 
the  results  of  the  last  two  Iterations. 

With  the  eigen  angle,  en,  for  mode  n  determined  the  following  mode  para¬ 
meters  are  readily  evaluated 


phase  velocity  at  the  ground  * 
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excitation  factor  as  defined  by  Wait  = 


(III-8) 


•  lit/  .  ~  ( I  +  >R\\)  (l  p2 

Mode  mixing  parameter  (printed  out  as  polarization)  = 


where 
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The  subscript  r  in  Eq.  (6)  implies  the  real  part  while  the  subscript  i  in 

Eq.  (7)  implies  the  imaginary  part.  The  subscript  on  q  represents  the  value 

of  z  at  which  q  is  evaluated.  For  example,  q*  means  that  Eq.  (15)  is  to  be 

2 

evaluated  for  Z  =0.  Similarly,  the  subscript  in  n  represents  the  value  of  Z 
for  which  Eq.  (4)  is  to  be  evaluated.  The  functions  h^  and  h^  are  the  modified 
Hankel  functions  of  order  1/3.  The  primes  in  Eqs.  (12)  and  (13)  denote  de¬ 
rivatives  with  respect  to  the  argument.  The  quantity  h’  which  occurs  in  Eq.  (8) 
is  called  REFL  HT  in  the  program  and  has  been  introduced  solely  for  the  purpose 

of  normalizing  the  excitation  factor  for  the  vertical  dipole  in  a  way  which  is 

(ll) 

consistent  with  Wait's  definition  '  .  Closely  related  quantities  also  in¬ 

cluded  in  the  print  out  are  termed  EXTRA  MAG  and  EXTRA  ANGLE.  These  quantities 
are  simply  the  mamitude  and  argument  of  the  right  hand  side  of  Eq.  (8)  with 
the  factor  (-i  ■  )  divided  out.  The  quantities  EXTRA  MAG  and  EXTRA  ANGLE  are 

used  in  conjunction  with  the  mode  sum  analysis  discussed  in  section  IV.  The 
quantity  printed  out  as  polarization,  that  is  the  right  hand  side  of  Eq.  (9), 
is  in  fact  not  a  measure  of  the  polarization  but  as  shown  in  reference  14  is 
very  closely  related  to  the  ratio  of  the  vertical  to  horizontal  components  of 
polarization  which  comprise  each  mode  under  the  general  condition  of  an  aniso¬ 
tropic  ionosphere.  The  right  hand  side  of  Eq.  (9)  becomes  meaningless  under 
conditions  of  an  isotropic  or  near  isotropic  ionosphere.  For  then  the  numerator 
of  equation  (9)  is  very  small  because  of  a  small  conversion  coefficient.  On 
the  other  hand  the  denominator  of  Eq.  (9)  even  for  a  vertically  polarized  mode 
will  not  vanish  in  general  because  of  round  off  error  so  that  the  entire  ratio 
is  then  dependent  upon  round  off  error  and  in  fact  can  have  magnitude  less  than 
unity  even  through  the  mode  is  vertically  polarized. 

Finally,  the  quantity  printed  out  under  the  label  THETA  PRIME,  Is  simply 
the  value  of  the  eigen  angle  at  the  ground  (recall  that  the  eigen  angle  which 
results  from  the  process  of  interation  is  referred  to  height  K). 
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IV.  CALCULATION  OF  RADIO  FIELD  STRENGTH 

For  the  frequencies  considered  for  waveguide  propagation  the  transmitting 
antenna  is  generally  a  vertical  monopole  above  the  surface  of  the  earth.  Be* 
cause  the  earth  is  a  good  conductor  the  radiation  characteristic  of  the  antenna 
is  virtually  that  of  a  dipole  radiating  into  a  half  space.  The  power  radiated 
Pr  is  then 


watts 


(IV-1) 


where  I  is  the  rms  current,  1  is  the  effective  height  of  the  antenna  and  X  *8 
the  radio  wavelength.  The  electric  field  at  a  distance  D  from  the  transmitter 
for  the  n  th  mode  can  be  written  as: 


sm 
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where  is  377  ohms,  c  is  the  speed  of  light,  f  is  radio  frequency,  a  is  the 
radius  of  the  earth,  Cn  is  the  cosine  of  0n  the  n  th  eigen  angle,  k  is  wave 
number,  and  x  in  the  magnitude  of  the  exitation  factor  (extra  mag.).  A 
practical  form  to  use  is 


where 


_C  oc 

Ziv/m/K ui)  =  C  e  X 


.2041 


C,  -  A 

C*=  -II52D 


(fsin&j 


(iv-3) 


and  a  *  attenuation  coefficient  in  db/me  game  ter. 

The  relative  phase  of  the  signal  (relative  with  respect  to  free  space  pro¬ 
pagation  or  to  a  local  oscillator)  is  given  by 
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(iv-M 


where  <p  is  the  relative  phase  at  the  receiver,  c  is  the  speed  of  light,  IT  is 
the  phase  velocity  and  cp*  is  the  argument  of  the  excitation  factor  in  radians. 
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V.  RUNNING  THE  PROGRAM 

A.  The  profile  deck(s) 

The  program,  whether  run  for  WAVEGUIDE  SOLUTIONS  or  for  REFLECTION 

COEFFICIENTS  can  be  run  with  or  without  ion  density  profiles.  The  profiles 

are  set  up  the  same  way  for  either  case.  The  electron  density  profile  comes 

first  and  is  preceded  by  a  SPECIES  1  PROFILE  card  and  an  alpha-numeric  card 

punched  with  information  needed  to  identify  the  run.  The  profile  starts  at 

the  top  of  the  ionosphere  and  contains  the  height  in  kilometers  and  the 

electron  density/cubic  centimeter  for  that  height.  The  height  is  punched  in 

columns  2-7  and  the  density  is  punched  in  columns  14-21.  The  decimals  are  in 

columns  5  and  15.  If  an  electron  density  of  zero  is  used  the  program  will 
— I4O 

change  it  to  10  .  Although  the  program  is  often  run  with  an  electron  den¬ 
sity  for  every  half  kilometer,  it  may  be  run  using  only  the  significant 

heights.  For  example,  only  two  cards  are  needed  for  an  exponential  profile: 

One  at  the  top  and  one  at  the  bottom  of  the  profile.  The  program  interpolates 
in  a  straight  line  between  the  logs  of  the  densities  and  heights  given  in  the 
profile  deck.  The  profile  deck  ends  with  a  card  punched  with  nines  in  the 
first  8  columns.  If  the  program  is  to  be  run  without  ions,  the  input  para¬ 
meters  immediately  follow  this  nines  card.  If  the  program  is  to  be  run  with 
ions,  this  nines  card  is  followed  by  a  SPECIES  2  PROFILE  card,  another  card 
with  identification,  then  the  top  of  the  positive  ion  density  profile  and  so 
on.  This  also  ends  with  a  nines  card  followed  by  a  SPECIES  3  PROFILE  card, 
an  identification  card,  and  the  negative  ion  density  profile.  Again  a  nines, 
and  then  the  input  parameters.  However,  since  the  negative  ion  density  is 

equal  to  the  positive  ion  density  less  the  electron  density  for  each  height, 

the  program  will  calculate  the  negative  ion  density  if  a  NEGIONS  card  is 
placed  at  the  end  of  the  positive  ion  density  profile.  The  ion  density  pro¬ 
files  must  have  the  same  heights  as  the  electron  density  profile.  A  sample 
profile  with  ions  included  is  shown  in  Table  A. 

B.  Waveguide  Solution  Input  Parameters 

The  input  parameters  are  placed  anywher?  on  the  card  but  according  to  the 
following  format.  The  names,  e.g.,  MAGFIELD,  must  have  no  spaces  between  the 
letters.  The  name  is  followed  by  one  space,  an  -  sign,  another  space,  and  then 
the  value  of  the  parameter.  If  a  parameter  has  more  than  one  value  they  are 
listed  on  the  card  or  on  following  cards  with  one  or  more  spaces  between  them. 
The  parameter  name  is  only  used  once. 
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Table  A 


The  approximate  method  will  find  up  to  25  solutions  in  10  minutes  and 
the  exact  method  might  find  one  solution  in  that  amount  of  time.  Therefore 
it  expediter  the  search  for  mode  solutions.  Often  these  solutions  are 
sufficiently  accurate  for  the  purpose.  If  this  is  not  the  case  they  can 
serve  as  excellent  starting  solutions  when  running  the  exact  method.  To 
set  up  an  approximate  run,  only  two  cards  must  be  added  to  an  exact  waveguide 
run.  These  are  the  RPOLY  card  and  the  TLIST  card  which  are  explained  below. 

AZIMUTH  is  the  direction  of  propagation  in  degrees,  measured  east  of  mag¬ 
netic  north.  West-to-east  propagation  corresponds  to  an  azimuth  of  90*  and 
east-to-west  propagation  corresponds  to  an  azimuth  of  270*. 

DIPANGLE  la  the  magnetic  dip  angle  in  degrees,  measured  from  the  vertical 
and  is  always  between  0  and  90  degrees  for  propagation  in  the  northern  himis- 
phere  and  between  90  and  180  degrees  for  propagation  in  the  southern  hemisphere. 

MAGFIELD  is  the  intensity  of  the  magnetic  field  in  Webers/square  meter. 

Use  0.0  for  no  mag  field. 

COEFFNU  and  EXFNU  are  the  values  used  to  find  ^  ,  the  collision  frequency 
where 


T)  =  COEFFNU  x  exp  (EXPNU  x  Z) 

and  Z  is  in  meters,  EXPNU  is  in  inverse  meters  and  COEFFNU  is  in  collisions/ 
second.  If  ions  are  Included,  a  separate  collision  frequency  is  needed  for 
each  species  and  is  listed  on  each  card  in  the  order  of  the  profiles.  For 
example  COEFFNU  =  4. 3O3EII  I.O76EIO  I.O76EIO  would  mean  that  the  coefficient 
of  the  collision  frequency  is  4.3O3EII  (U.  3O3  x  1011)  for  the  electron  density 
profile  and  I.076EIO  for  the  two  ion  density  profiles. 

EPSILON  and  EPSILONO  are  the  permittivities  of  the  ground  and  free 
space  respectively  and  are  to  be  expressed  in  units  of  farads/meter. 

D  is  any  height  below  which  ionospheric  effects  can  be  considered  neglig¬ 
ible  and  must  be  consistent  with  the  bottom  of  the  profile  deck.  Care  must 
be  taken  in  satisfying  the  preceding  condition.  If  there  is  doubt  as  to 
whether  a  suitable  value  has  been  chosen,  it  is  advisable  to  run  the  program 
to  different  heights,  D,  to  verify  that  the  eigen-angles  have  stabilized. 
Generally,  stability  to  better  than  +  0.05°  in  the  ground  eigenangle  is  ob¬ 
tained. 
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H  is  the  height,  expressed  in  kilometers,  which  represents  physically 
the  height  at  which  the  modified  refractive  index  becomes  unity  and  is  the 
height  at  which  the  unprimed  eigenangle  is  measured.  Theory,  good  to  first 
order  in  the  parameters  h/a  and  z/a,  predicts  that  the  final  ground  solutions 
are  independent  of  H.  However,  investigation  has  shown  that  the  results  vary 
somewhat  when  H  changes  from  90  km  to  the  ground  (by  approximately  a  tenth  of 
a  degree  in  the  ground  eigen  angle).  This  variation  is  attributed  to  higher 
order  effects  in  the  earth  curvature  correction.  Although  contributing  to  some 
uncertainty  in  the  final  results,  these  variations  can  be  considered  small  and 
are  sufficiently  accurate  for  most  purposes.  Generally,  it  is  sufficient  to 
choose  both  D  and  H  to  be  the  height  on  the  last  card  of  the  electron  density 
profile. 

ALPHA  is  in  inverse  kilometers  and  is  equal  to  2/(earth*s  radius  in  kilo¬ 
meters).  If  flat  earth  programs  are  run,  a  very  small  number  (e.g.  10”^  km-1) 
should  be  used  instead  of  zero. 

RELPREC  is  a  parameter  which  helps  determine  how  accurately  the  calculations 
are  to  be  done.  A  RELPREC  of  1.0  would  be  fast  but  not  very  accurate,  2.0  is 
usual,  and  3.0  or  4.0  would  give  about  maximum  accuracy. 

SIGMA  is  the  ground  conductivity  expressed  in  units  of  mhos  /me  ter. 

REFLHT  is  chosen  to  be  a  reasonable  estimate  of  where  reflection  takes 
place.  It  is  used  in  both  the  approximate  and  the  exact  methods  as  a  factor 
in  obtaining  WAITS  EXCITATION  FACTOR.  It  is  also  used  in  obtaining  the  approxi¬ 
mate  solutions.  It  is  in  kilometers. 

FREQKHZ  is  the  radio  frequency  in  kilohertz/second. 

THETAINC  is  a  card  used  to  control  the  iterations.  A  THETAINC  of  1.0 
insures  that  each  succeeding  Iteration  will  be  no  more  than  one  degree  (in 
real  and  in  imaginary  parts)  from  the  preceeding  one.  This  keeps  the  program 
from  arriving  at  a  solution  other  than  the  one  desired.  However,  if  little 
is  known  about  the  eigen  angle  solutions,  a  THETAINC  of  5.0  or  even  greater 
would  be  best. 

MRATIO  is  used  only  when  ions  are  included.  The  first  number  Is  unity; 
the  second  the  ratio  of  positive  ion  mass  to  electron  mass;  and  the  third,  the 
ratio  of  negative  ion  mass  to  electron  mass. 

CHARGE  is  also  necessary  only  when  ions  are  included.  The  first  number  is 
the  charge  of  electrons,  -1.0,  the  second,  the  charge  of  positive  ions,  +  1.0, 
and  the  third  the  charge  of  negative  ions,  -1.0. 
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RPOLY  =  1  if  approximate  solutions  are  desired. 

TLIST  is  a  list  of  complex  angles  used  in  the  approximate  method  to 
obtain  exact  reflection  coefficients  which  are  used  in  an  interpolation 
scheme  to  find  the  eigenangles  of  the  problem.  The  closer  the  elements  of 
the  TUST  are  to  the  final  eigen  angle  solutions  the  closer  the  solutions 
will  be  to  exact  solutions.  For  example,  if  all  modes  were  derived  between 
70*  and  80“,  the  TLIST  might  look  like  the  following: 

TUST  =  70  -3  75  -2  80  -1 

Then  all  eigens  in  that  range  could  be  run  in  very  little  computer  time. 

E.g.  the  EIGEN  card  might  appear  as 
EIGEN  =  70.6  -2.7  73.5  -1.2  7k  -1  78  -2.5 

If  that  is  not  accurate  enough,  and  an  approximate  solution  near  (76,  -1.6)*  is 
desired  the  TUST  might  look  like  the  following: 

TUST  =  75  -1.8  76  -1.6  77  -l.U 

EIGEN  precedes  the  trial  eigen  angles.  As  many  eigen  angles  as  desired 
are  punched  on  cards  in  the  format  of  real  part,  imaginary  part,  real  pert, 
imaginary  part,  and  so  on.  Each  number  separated  by  one  or  more  spaces.  If 
more  than  one  card  is  needed  the  numbers  are  Just  continued  on  the  following 
cards  without  reusing  the  same  EIGEN. 

WAV2GUID  is  the  command  to  the  computer  to  start  and  complete  the  wave¬ 
guide  solutions. 

See  Table  B  for  a  sample  approximate  waveguide  run  with  ions  included. 

These  parameters  would  follow  the  profile  decks.  For  the  run  to  be  the 
exact  method,  only  the  RPOLY  card  and  the  TUST  card  need  be  removed. 

C.  Reflection  Coefficient  Input  Parameters 

FREQKHZ  (See  Waveguide  input  parameters) 

AZIMUTH 
DIPANG LE 
MAGFIELD 
COEFFNU 
EXFNU 
D 
H 

ALPHA 
RELPREC 


II 

•I 

II 

II 

!• 

If 

II 

II 

ft 
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Table  B 


A7Imuth  a  8  3.0 

D!PA*r,l_E  =  17  C 
M  A  R  f I P  L  U  =  5,0E-(ic 

COgrr'  ,.i  a  4.3l'3Fll  1.C76F10  1.&74F10 

EVPMU  a  -1.622E-04  -1.622E-0 4  -1.62?E-04 

EPSTLPN  =  7.17201^4F-10 
PPSTLONO  =  B.B5434F-1? 

D  a  P.n 
H  a  C1.0 

AI.PWa  s  3  ,  i 4 1 -  C 4 
PF|_P*FC  a  4.0 
STRf'A  a  4.64 
RPFI  HT  a  ftp,  ,  c 

MP  a  t  i  n  =  i,o  SBCTo.C  *’6000.0 
ruAPRP  -  “i,c  i .«•  -l.o 
TMETaTWC  a  l.C 
RPfll  V  a  1 

tl  ICT  a  7^.?  -7.0  8? .  0  -k’.O  80.  *5  -1.0 

FT3PN  a  75, i-  -2.0  76.0  -2.0  77.0  -2.0  78.0  -2-0  79.0  ->.0  80.0  -2.0 
8?.P  -2.0  83.0  -2.0  fc4,o  -2.C  85.0  -2.0  86.0  -2.0  87.0  -2.0  84.0  -? 
87 . *  -2.C  &H.5  -2.0  6P.5  -2.0 
fppr Z  *  lf.0 
WAVM?l'!D 
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MRATIO  (if  ions  are  in¬ 
cluded) 

CHARGE  "  '• 

THETA  is  the  complex  angle  of  incidence  in  degrees. 

See  Table  C  for  a  sample  input  for  reflection  coefficients  for  electrons 
only.  The  reflection  coefficient  program  only  runs  exact. 

RC  is  the  command  to  the  computer  to  start  and  complete  the  reflection 
coefficients. 

D.  Linking  Several  Computer  Runs 

By  doing  several  different  runs  in  one  computer  input  one  can  sometimes 
save  on  turnaround  time  and  cut  down  on  card  duplication  of  either  inputs  or 
profiles. 

More  than  one  data  deck  is  run  at  a  time  by  following  the  WAVEGUID.  card 
with  a  SPECIES  1  PROFILE  card  and  so  on  through  the  profile  or  through  the 
ion  profiles  as  before.  Doing  this  clobbers  the  profiles  which  were  originally 
in  the  computer  but  all  the  other  input  remains  as  it  was.  At  the  next  WAVE¬ 
GUID  card  the  program  starts  on  this  "second"  program.  If  any  changes  In  input 
parameters  were  desired  on  this  second  run,  they  would  be  changed  between  the 
profile  and  the  WAVEGUID  card.  For  example,  if  D  =  0.0  on  the  first  profile 
and  D  =  65.O  on  the  second,  the  new  profile  would  end  with  a  nines  card  and 
then  NEGIONS,  if  ions  were  included,  then  D  =  65. 0  and  then  WAVEGUID. 

Changes  in  the  parameters  for  the  data  deck  are  stated  between  WAVEGUID 
cards.  All  other  inputs  would  remain  unchanged.  For  example,  if  several  mag 
fields  are  desired  for  a  given  profile  and  set  of  parameters,  only  mag  field 
would  be  changed  before  calling  WAVEGUID  again.  For  example 

MAGFIELD  =  U.0E-05 
WAVEGUID 
MAGFIELD  =0.0 
WAVEGUID 

It  is  sometimes  justifiable  under  disturbed  conditions  such  as  SID's  or 
man-made  ionization  anomalies  to  neglect  the  influence  of  the  earth's  magnetic 
field.  The  modal  equation  then  factors  into  a  product,  one  factor  of  which 
contains  solutions  for  vertical  polarization,  while  the  other  factor  contains 
solutions  for  horizontal  polarization.  A  provision  is  thus  made  whereby  under 
those  circumstances  either  the  vertical  or  the  horizontal  solutions  may  be 
obtained  independently.  Of  course  for  vertical  dipole  excitation,  the  verti¬ 
cally  polarized  modes  ar*  the  physically  meaningful  modes.  So  when  vertical 
modes  are  desired  and  the  profile  is  run  with  no  mag  field,  the  command  TYPEITEP=  1 
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Table  C 


SDE'1IPS  1  profile 

sample  PROFILE  FLPCTRONf5  0 . J L v 


tOP.no 

7.02F  04 

?*'.  CTO 

3 . 0  6  p  0  4 

90 .0(1 

1.33E  04 

05.0  0 

6.12P  03 

on  .no 

2.80F  03 

7*.  oo 

1.34E  03 

70  .  no 

5.72M  02 

65.0  0' 

2.07 F  02 

60.0  0 

3 , 4  0  c  01 

55.00 

6.50E  00 

50 . 00 

4 .37F-01 

45,00 

1 . 6  0  E  -  0 1 

4  0 . 0  0 

5.42F-02 

35.07V 

1.58P-02 

30.00 

4 . 03E-03 

25.no 

9 , 54p:_Q4 

20 . 00 

2. 19P-04 

15.00 

5.32E-05 

in.no 

1 . 4  7  p  -  0  5 

5. 00 

5 . R3F-06 

0 .00 

O.OPF  00 

99909099 

A7  I  M(jj  h  =  33. 

a 

DTPA'MPLE  =  17 

.  J 

MAif i eld  =  5. 

0  E  -  0  5 

EYPX'J  s- 1.6??=- 04 
COeFFN'U  =  4.303E11 
D  =  0.0 
H  =  0.0 

ALPMA  s  3.14E-04 
RFLPREC  =  2.0 
rpgOKF?  =  15 ;o 
T  MET  A  =  30.5  -1.5 
RC 


20 


will  force  only  parallel  solutions.  If  desired,  horizontal  solutions  only 
will  result  from  the  command  TYPEITER  =  2.  Either  command  may  be  used  when 
propagation  E-W  or  W-E  is  run  at  the  geomagnetic  equator  (i.e.  propagation 
direction  transverse  to  the  magnetic  field). 

E.  Waveguide  Solution  Computer  Printout 

Line  Explanation  of  Table  D 

1,  2  Real  and  imaginary  parts  of  the  eigen  angle  followed  by  their 

reflection  coefficients. 

3  Difference  between  final  two  eigen  angles. 

4  Real  and  imaginary  parts  of  the  F  function.  See  III. 

5  Phase  velocity  of  mode  at  the  earth's  surface. 

6  Phase  velocity  divided  by  the  speed  of  light. 

7  Attenuation  of  the  vertical  E-field  in  dB/1000  km  at  the  ground. 

3,  9  Magnitude  and  angle  of  Wait's  excitation  factor  which  has  been 

normalized  by  dividing  out  the  -ikh'/2  term. 

10,  11  Wait's  excitation  factor  explained  in  (8).  These  excitation  factors 
are  for  a  vertical  dipole  antenna.  See  III -8. 

12,  13  These  give  an  indication  if  the  eigen  angle  is  horizontally  or 
vertically  polarized  and  are  explained  in  III-9. 

14  The  eigen  angle  at  the  ground. 

F.  Reflection  Coefficient  Computer  Printout 
Explanation  of  Table  E. 

The  computer  prints  out  the  height  first  then  across  the  page,  the  modulus 
and  argument  in  radians  of  „ R,,  ,  XRA  ,  XR„  ,  (|RX  .  At  the  bottom  of  the 
printout  it  prints  XR„  /  XRX  and  x0„  -  j9x  . 
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APPENDIX  A 


Some  Notes  on  Organization  of  the  Waveguide  Program. 

a.  Functional  Dependence  Chart 

The  general  features  of  the  waveguide  program  are  illustrated  on  the  next 
page  in  a  kind  of  "functional  dependence"  chart  (it  is  not  a  flow  chart). 
Control  in  the  program  begins  in  PROGRAM  INRJT  which  calls  on  the  waveguide 
control  routine,  which  in  turn  calls  on  the  iteration  routine,  etc.  Each 
routine  is  initiated  by  the  one  above  it  on  the  chart.  The  items  in  the 
chart  correspond  roughly  to  Fortran  subprograms. 

b.  Program  Input  and  Dummy 

Program  input  is  used  at  NEDC  for  various  program  combinations  other  than 
the  waveguide  program  described  in  this  report.  Hence  some  of  the  programming 
in  this  routine  is  not  applicable  to  the  waveguide  program.  Names  of  routines 
called  on  by  INRJT,  but  not  in  this  program  deck,  are  programmed  as  entry  points 
in  DUMMY  to  make  compilation  possible. 

The  omitted  programming  consists  of  an  alternate  full-wave  solution  (a 
modified  Inoue -Horowitz  Integration)  and  a  scheme  for  fitting  electron  density 
profiles  to  reflection  coefficient  data.  Description  of  these  routines  is 
beyond  the  scope  of  this  report. 

The  program  calls  on  ZERO  CLK  and  PRINTIME.  These  are  written  in  CDC 
106k  assembly  code  and  are  used  to  time  various  segments  of  the  program.  It 
is  probably  best  that  the  user  substitute  his  own  clock  routine. 


2k 


t 


PROGRAM  INPUT 


COMMON/INPUT/THETA, FREQ  KHZ » AZ I MUTH.DI P  ANGLE  *MAG  FIElD*ALPHA,H* 

S  O.REL  PREC »TOP  HT, LOWEST  HT.IN  EXTRA(88) 

COMMON/SP  INPUT/NR  SPEC*COEFF  NU(3)»EXP  NU( 3 ) *  CHARGE! 3) *M  RATIO(3) 
COMMON/FLG  INPUT/DEBUG. RC  PRINT»POLY  FLAG*INTEG  UP»WF  FLAG* 

S  NEG  ELECT, FLG  EXTRA(19) 

COMMON/WG  INPUT/TYPE  ITER, MAX  DELTA, EPSILON, EPSILON  0, SIGMA, 

$  WGIN  SKIP, 

$  REEL  HT  *R  POLY,T  L I S T<1 1 ) , E I  GEN ( 30 ) • THETA  INC, 

$  T,PHI ,DTHETA,WG  EXTRA ( 5 ) 

COMMON/PF  INPUT/NR  V»NR  F, 

$  EXPER  F(6)*£XPER  VAL < 24 ), EXPER  UNC(24), 

$  DA, RMS  TOL ,NR  TIMES, NR  FRACT , FRACT ( 7 ) ,DA  FACTOR, 

$  PF  EXTRA ( 32 ) 

COMMON/PR  INPUT/NR  A,HTO»HT  UNIT»A(9)*NR  PTS.FN  FACTOR, 
s  pr  extra ( i i j 

COMMON/l/SKIPl ,HT  L I  ST ( 2 56 ) ,LOG  N  LIST(256*3) 

NAMELIST/DATA/THETA. FREQ  KHZ, AZIMUTH, DIP  ANGLE* MAG  FIELD, ALPHA, H, 


$  D,REL  PREC, TOP  HT, LOWEST  HT , 

$  COEFF  NU.EXP  NU, CHARGE, M  RAT  1 0 .SPEC  I ES , 

S  DEBUG, POLY  FLAG, NEG  ELECT, 

S  TYPE  ITER, MAX  DELTA, EPSILON, EPSILON  0, SIGMA, 

S  REFL  HT,R  POLY,T  LIST, EIGEN, THETA  INC, T, PHI, 

S  NR  V *NR  F, EXPER  F, EXPER  VAL, EXPER  UNC, 

S  DA, RMS  TOL, NR  TIMES, NR  FRACT , FRACT , DA  FACTOR, 

$  NR  A,HTO»HT  UNIT  ,A ,FN  FACTOR 


COMPLEX  THETA 

COMPLEX  EIGEN, T  LIST  »DTHET  A 
REAL  MAG  FIELD, LOWEST  HT ,M  RATIO 
REAL  MAX  DELTA 
REAL  LOG  N  LIST 

INTEGER  DEBUG, RC  PRINT, POLY  FLAG.WF  FLAG 
INTEGER  ALPHANUM, SPECIES 
INTEGER  TYPE  ITER.R  POLY 
DIMENSION  ALPHANUM(IO) 


DATA(ALPHA«0.0) ,(H*0.0) 

DATA  ( D=0,0 ) 

DATA ( REL  PREC-3,5) 

DATA (LOWEST  HT«0.0) 

DAT  A ( SPEC  I ES* 1 ) , ( NR  SPEC«1) 

DATA ( M  RAT  10(1 1 ■ 1*0) .(CHARGE ( 1 )«-1.0) 
DATA(DEBUG»0) 

DATA ( RC  PRINT«1) 

DATA (POLY  FLAG=0) 

DATA ( INTEG  UP*0) 

DATA(WF  FLAG«0) 

DATA ( NEG  ELECTRO) 

DATA ( TYPE  ITER*0 ) , (MAX  DELTA«0.0) 

DATA ( R  POLY*0 ) 

DATA(  T  LIST«11(  (O.O.O.OM  ) 

DATA (EIGEN«30(( 0.0, 0.0) )  ) 

DATAtTHETA  INC*0.1) 
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DATA (DTHETA=< 0.05 *0.01 )  ) 


DATAtNR  V  =  9 ) 

DATA ( NR  F  =  3) 

DATA(DA=1.0) 

DATAIRMS  TOL  =  1 • 0  5  ) 

DAT  A ( NR  T I MES= 1 ) 

DATA (NR  FRACT*7 ) 

DATA (FRACT=0.015625 *0.03125 *0.0625 *0.125 *0.25 *0.5* 1.0) 
DATA ( DA  FACTOR3 1.0) 

DAT  A ( NR  A  =  6) 

DATA ( HTO=60.0)*(HT  UNIT*15.0) 

DAT  A ( A  *9 ( 0 .0 ) ) 

DATA ( FN  FACTOR=1.0E03) 


PRINT  100 
100  FORMAT ( 1H1 ) 

10  NW  =  INLIST(DATA*N) 


IF  ( NW 

•  EO. 

7HPR0FILE) 

GO  TO 

20 

IFINW 

•  EO. 

7HREVERSE) 

GO  TO 

30 

IF  ( NW 

•  EO. 

7HNEGI0NS) 

GO  TO 

35 

IF  ( NW 

•  EQ. 

2HRC)  GO  TO 

40 

IF  ( NW 

•  EQ. 

8HINITIALA) 

GO  TO 

45 

I F  ( NW 

•  EQ. 

2HWF)  GO  TO 

i  50 

IF  ( NW 

•  EO. 

8HWAVEGUID) 

GO  TO 

60 

I F  ( NW 

•  EQ. 

6HENC0EF)  GO  TO  70 

I F  ( NW 

•  EO. 

7HFITPR0F) 

GO  TO 

71 

I F  ( NW 

•  EO. 

8HENUNCERT) 

GO  TO 

72 

IFINW 

•  EO. 

4HQUIT)  GO 

TO  99 

GO  TO 

90 

20  READ  201*(ALPHANUM( J) *J«1 *10) 

PRINT  202  •  ( ALPHAN'JM(  J)  *J=1*10) 

K  =  SPECIES 

I F ( K  .GT.  NR  SPEC)  NR  SPEC  *  K 
L  =  1 

21  READ  201 *(ALPHANUM( J) *J*1 .10) 

?01  FORMAT ( 10A8 ) 

PRINT  202  * ( ALPHANUM( J)*J=1*10) 

202  FORMAT ( 1H  »10A8) 

IF(ALPHANUM( 1)  .EO.  8H99999999)  GO  TO  22 
DECODE (  21  *  203 .ALPHANUM )  HT  *EN 

203  FORMAT (F7.2.5X.E9.2) 

IF  I K  .NE.  1  .AND.  HT  .NE.  HT  LIST(L))  GO  TO  90 
HT  LIST(L)  3  HT 

I F  ( L  .NE.  1  .AND.  HT  LIST(L)  .GE.  HT  LIST(L-D)  GO  TO  90 
IF ( EN  .EO.  0.0)  EN  =  1.0E-40 
LOG  N  L I  ST  I L  *K )  «  LOGF(EN) 

L  =  L  +  l 
GO  TO  21 

22  TOP  HT  =  HT  LIST(l) 

LOWEST  HT  3  HT  LlST(L-l) 

NR  PTS  *  L-l 

INTEG  UP  *  0 
GO  TO  10 
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30  L  HALF  =  NR  PTS/2 
DO  31  L=1 »L  HALF 
M  =  NR  PTS+l-L 

TEMP  =  HT  LIST(L)  S  HT  LIST(L>  «  HT  LIST(M)  S  HT  LlST(M)  =  TEMP 
TEMP  =  LOG  N  L I  ST ( L* 1 )  S  LOG  N  LISTIL.1)  =  LOG  N  LIST(M*1> 

LOG  N  L 1ST (M.l )  *  TEMP 

31  CONTINUE 

DO  32  L*1 *NR  PTS 

32  HT  LISTCL)  =  -HT  LIST(L) 

INTEG  UP  =  1-INTEG  UP 

TEMP  =  TOP  HT  S  TOP  HT  «  -LOWEST  HT  S  LOWEST  HT  ■  -TEMP 
GO  TO  10 

35  DO  36  L=1 *NR  PTS 

EN  =  EXPF  (  LOG  N  L I  ST  ( t  .2  )) -EXPF  (  LOG  N  LIST(L*1)> 

IF ( EN  .EQ.  0.0)  EN  =  1.0E-40 
LOG  N  L I  ST ( L  »  3 )  =  LOGF(EN) 

36  PRINT  301 *HT  LIST(L)*EN 
301  FORMAT ( 1H  *F7. 2 ♦ 5X .E9.2 ) 

NR  SPEC  *  3 
GO  TO  10 

C  RC 

40  CALL  ZERO  CLK 
D  »  LOWEST  HT 
CALL  INTEG 
CALL  PR  I  NT  I  ME 
PRINT  100 
GO  TO  10 

45  D  «  LOWEST  HT 
CALL  INITIAL  A 
GO  TO  10 

C  WAVEFIELDS 

50  CALL  ZERO  CLK 
WF  FLAG  «  1 
D  «  LOWEST  HT 
CALL  INTEG 
CALL  BACKWARD 
WF  FLAG  «  0 
CALL  PRINTIME 
GO  TO  10 

C  WAVEGUID 

60  CALL  WAVEGUID 
GO  TO  10 

C  PROFILE  FITTING 
TO  CALL  EN  COEF 
GO  TO  10 

71  D  *  0.0 

CALL  FIT  PROF 
GO  TO  10 

72  CALL  EN  UNCERT 
GO  TO  10 
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C  ERROR  EXIT 
90  PRINT  900 

900  FORMAT ( 1H0 ♦ 18HERR0R  IN  DATA  DECK) 
99  CONTINUE 
END 


SUBROUTINE  DUMMY 
C  FOR  BUDDEN  WAVEGUIDE 

ENTRY  FIT  PROF 
ENTRY  EN  UNCERT 
ENTRY  EN  COEF 
ENTRY  BACKWARD 
ENTRY  INITIAL  A 
RETURN 

END 
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FUNCTION  CPLX  FCT(ARG) 


COMPLEX  CPLX  FC T ♦ ARG » I 
REAL  NEG 

DATA( SAVE  REAL* 1 .0E99 ) • ( SAVE  IMAG=1.0E99) 
DATA(  IMO.Oil.OII 
DATA (PI =3.1415926536) 


ENTRY  CPLX  SORT 

ARG  REAL  «  ARG  $  ARG  IMAG  =  -I*ARG 

RHO  =  SQRTF(ARG  REAL*ARG  REAL+ARG  IMAG*ARG  IMAG) 

ABS  REAL  =  ABSF ( ARG  REAL) 

I F  < RHO  .LT.  ABS  REAL)  RHO  =  ABS  REAL 
I F { ARG  IMAG  .GF.  O.O)  1.2 

1  CPLX  FCT  =  SORTF  (  (  RHO+ARG  REAL  )  *0  •  5  )  +  1  *SQRTF  (  <  RHO-ARG  REAL)*0.5) 
RETURN 

2  CPLX  FCT  *  -SQRTF( (RHO+ARG  REAL ) *0 • 5]+ I *SQRTF (( RHO-ARG  REAL)*0.5) 

RETURN  _  -  -  -  - -  ' 

ENTRY  CPLX  LOGF 

ARG  REAL  *  ARG  *  ARG  IMAG  =  - I *ARG 

RHO  =  SQRTF ( ARG  REAL*ARG  REAL+ARG  I MAG*ARG  IMAG) 

I F ( RHO  .EQ.  0.0)  GO  TO  13 
AeS  REAL  =  ABSF ( ARG  REAL) 

I F { RHO  .LT.  ABS  REAL)  RHO  =  ABS  REAL 
I F ( ARG  IMAG  .GF.  0.0)  11,12 

11  CPLX  FCT  =  LOGF (RHO)+I*ACOSF( ARG  REAL /RHO )  t  RETURN 

12  CPLX  FCT  *  LOGF < RHO )+I*(2.0*PI-ACOSF( ARG  REAL/RHO)) 

RETURN 

13  CPLX  FCT  =  -200.0  5  RETURN 
ENTRY  CPLX  EXPF 

ARG  REAL  «  ARG  *  ARG  IMAG  =  -I*ARG 
IF ( ARG  IMAG  .FO.  0.0)21,22 

21  CPLX  FCT  =  EXPF ( ARG  REAL)  *  RETURN 

22  IF ( ARG  REAL  .EQ.  0.0)23,24 

23  CPLX  FCT  =  COSF(ARG  IMAG ) +1 *S I NF ( ARG  IMAG)  S  RETURN 

24  CPLX  FCT  «  EXPF(ARG  RE AL ) * < COSF < ARG  I  MAG ) ♦ I *S I NF ( ARG  IMAG)) 

RETURN 

ENTRY  CPLX  COSF 

ARG  REAL  =  ARG  $  ARG  IMAG  =  - I *ARG 
IF ( ARG  IMAG  .EQ.  0.0)31,32 

31  CPLX  FCT  =  COSF ( ARG  REAL)  S  RETURN 

32  IF ( ARG  REAL  .EQ.  SAVE  REAL  .AND.  ARG  IMAG  .EQ.  SAVE  IMAG)  GO  TO  33 
SAVE  REAL  =  ARG  REAL  S  SAVE  IMAG  =  ARG  IMAG 

COS  =  COSF (ARG  REAL) 

SIN  =  SINF ( ARG  REAL) 

POS  *  EXPF ( ARG  IMAG)  $  NEG  =  l.O/POS 
COSH  =  ( POS+NEG ) »0. 5 
SINH  =  ( POS-NEG ) *0. 5 

33  CPLX  FCT  =  C0S*C0SH-I*SIN*SINH 
RETURN 

ENTRY  CPLX  SINF 

ARG  REAL  =  ARG  t  ARG  IMAG  *  - I *ARG 
IF ( ARG  IMAG  .EQ.  0.0)41,42 


■>1 


41  CPLX  FCT  =  SINF(ARG  REAL)  S  RETURN 

42  IF { ARG  REAL  .EO.  SAVE  REAL  • AMD*  ARG  IMAG  .EQ*  SAVE  IMAG)  GO  TO  43 
SAVE  REAL  *  ARG  REAL  S  SAVE  IMAG  «  ARG  IMAG 

COS  =  COSF ( ARG  REAL) 

SIN  =  SINF ( ARG  REAL) 

POS  =  EXPFlARG  IMAG)  S  NEG  *  l.O/POS 
COSH  =  ( POS+NEG ) *0*  5 
SINH  =  (POS-NEG)#0.5 

43  CPLX  FCT  =  SI N*COSH+ I *COS*S I NH 
RETURN 

ENTRY  CONJ 

ARG  REAL  *  ARG  S  ARG  IMAG  *  -I* ARG 
CPLX  FCT  *  ARG  REAL-I*ARG  IMAG 
RETURN 

END 
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SUBROUTINE  MAG  ANGLE ( ARG* MAG* ANGLE ) 


COMPLEX  ARG,  I 
REAL  MAG 

DATA( I=(0. 0,1.0) ) 

DATA( RAD  TO  DEG=57. 29577951 ) 


ARG  REAL  «  ARG  $  ARG  IMAG  =  - I *ARG 

MAG  =  SORTF(ARG  REAL*ARG  REAL+ARG  IMAG*ARG  IMAG) 

COS  *  ARG  REAL /MAG 

I F ( COS  .GT.  1.0  .AND.  COS  .LT.  1.01)  COS  =  1.0 
IF ( COS  .LT.  -1.0  .AND.  COS  .GT.  -1.01)  COS  =  -1.0 
ANGLE  =  ACOSF ( COS)*RAD  TO  DEG 
IF ( ARG  IMAG  .LT.  0.0)  ANGLE  *  360.0-ANGLE 
RETURN 

END 


SUBROUTINE  QUART IC( B3*B2*B1*B0  *0* DEBUG) 


COMPLEX  B3  *32  »B1 »BO*Q» 

S  B3  SQ*H* I *G*H  PRIME  *  G  PRIME* 

S  SO  ROOT  *P  PLUS  *P  *LOG  P. 

S  CUBE  RTO  *CUBE  RT1*CUBE  RT2  *0MEGA1 *0MEGA2 * 

$  ROOT  P.ROOT  Q*ROOT  R. 

S  FUNCT I  ON  *CPLX  I* 

S  CPLX  SORT *CPLX  LOGF*CPLX  EXPF 

INTEGER  DEBUG. ERR  FLAG 
REAL  MAG  P'.US.MAG  MI NUS.MAG  F 
DIMENSION  0 ( 4 ) *P  RI ( 2 ) *FUNCT ION ( 4) 

DATA (0MEGA1= (-0.5* 0*8660254038) I » (OMEGA2® (-0*5 *-0*8660254038) ) 
DATA  ( CPLX  IMO.OihOII 
DATA ( T0L*1 .0E-02 ) 

EOUlVALENCE(P*P  RI) 


B3  SO  *  B3**2 
H  *  B2-B3  SQ 

I  *  B0-4.0*B3*B1+3.0*B2**2 
G  *  B1+B3* (-3.0*82+2.0*83  SO) 

H  PRIME  *  -1/12*0 

G  PRIME  «  -G**2/4.0-H#(H**2+3.0*H  PRIME) 

SQ  ROOT  «  CPLX  SORT ( G  PR IME»*2+4.0*H  PRIME**3) 

P  «  ( -G  PRIME+SQ  ROOT )*0.5 

MAG  PLUS  *  ABSF  ( P  RI  ( 1 )  )  ♦ABSF  (  P  RI(2)) 

P  PLUS  *  P 

P  *  (-G  PRIME-SQ  ROOT )*0*5 

MAG  MINUS  «  ABSF ( P  R I ( 1 ) )+ABSF (P  RI ( 2  )  ) 

IF (MAG  PLUS  *GT.  MAG  MINUS)  P  «  P  PLUS 
LOG  P  *  CPLX  LOGF(P) 

CUBE  RTO  ■  CPLX  EXPFILOG  P/3.0) 

CUBE  RT1  =  0MEGA1*CUBE  RTO 
CUBE  RT2  *  0MEGA2*CUBE  RTO 

ROOT  P  ■  CPLX  SORT ( CUBE  RTO-H  PRIME/CUBE  RTO-H) 

ROOT  Q  «  CPLX  SORT (CUBE  RT1-H  PR  I  ME /CUBE  RT1-H) 

ROOT  R  ■  CPLX  SORT ( CUBE  RT2-H  PRIME/CUBE  RT2-H) 

SIGN  •  -ROOT  P*ROOT  Q*ROOT  R*2.0/G 
IF (SIGN  *LT*  0.0)  ROOT  R  «  -ROOT  R 
0(1)  *  ♦ROOT  P+ROOT  Q+ROOT  R-B3 
0(2)  *  ♦ROOT  P-ROOT  Q-ROOT  R-B3 
Q( 3 )  «  -ROOT  P+ROOT  Q-ROOT  R-B3 
0(4)  *  -ROOT  P-ROOT  ©♦ROOT  R-B3 

DO  32  N*1 *4 

QR  *  Q(N )  S  ABS  QR  «  ABSF(QR) 

01  «  -CPLX  I *0 ( N )  5  ABS  QI  *  ABSF(QI) 

IF ( ABS  QR*TOL  .LT.  ABS  QI  )  GO  TO  31 

QI  ■  CPLX  I*( ( ( (QR+4.0*B3)*QR+6.0*B2 )*QR^4.0*B1)*QR^B0) 

S  /( ( (4.0*QR^12.0*B3)*QR+12.0*B2)*QR+4.0*B1 ) 

31  IF ( ABS  QI*TOL  .LT.  ABS  QR )  GO  TO  32 

OR  =  ( ( ( ( QI -CPLX  I*4.0*B3)*QI-6.0*B2)*QI+CPLX  I*4.0*B1 )*QI+BO ) 

S  /(((CPLX  I *4.0*01 +12 .0*83 ) *01 -CPLX  1*12 .0*B2 )*QI-4.0*B1 ) 

32  Q(N)  =  OR^CPLX  1*01 


IF ( DEBUG  #EQ.  0)  RETURN 


ERR  FLAG  *  0 
DO  51  N*1 #4 

IF ( DEBUG  .GE.  2)  PRINT  50l«N«0(N) 

501  FORMAT  (1H+.5HR00T  tU<3H  «  .C  (  El  2 .5  *  E13 .5  )  ) 

FUNCT  ION ( N )  =  C  (  <  < Q( N )*4.0*B3 > *Q ( N M-6.0*B2 )*Q( N > +4.0*B1 > *Q< N) +B0 ) 
f  /BO 

IF ( DEBUG  .GE.  2)  PRINT  502.N* FUNCT ION( N) 

502  FORMAT ( 1H  *40X * 9HFUNCT ION  •  1 1  * C ( E 1 1 • 3 »E1 1 • 3 ) > 

CALL  MAG  ANGLE! FUNCT I ONI N) *MAG  F. ANGLE  F) 

IF (MAG  F  .LT.  0.10)  GO  TO  51 
PRINT  503  *N#MAG  F 

503  FORMAT ( 1H  .5HR00T  »I1*23H  IS  NO  GOOD  MAG  F  »  .E10.3) 

ERR  FLAG  *  1 

51  CONTINUE 

IF ( ERR  FLAG  .EG.  1)  STOP 
RETURN 

END 
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SUBROUTINE  WAVEGUID 


COMMON /INPUT /THETA* IN  OMIT (98) 

COMMON /F LG  INPUT/FLG  SKIP.RC  PRINT *FLG  0MIT(23) 

COMMON/WG  INPUT/WGIN  SKIP(7)*R  POLY.WGIN  OMI T < 22  )  *EIGEN< 30) • 
S  WGIN  N0(10) 

COMPLEX  THETA*EIGEN 
INTEGER  RC  PRINT *R  POLY 


RC  PRINT  «  0 

IF ( R  POLY  .EG.  1)  CALL  GEN  R  POLY 
INDEX  «  1 

11  THETA  *  EIGEN( INDEX) 

CALL  ITERATE 

CALL  COMP  PROC 

EIGEN  RL  *  EIGEN( IN0EX+1 ) 

IF ( EIGEN  RL  .EQ.  0*0)  GO  TO  21 
INDEX  «  INDEX+1 
GO  TO  11 

21  RC  PRINT  «  1 
RETURN 

END 
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SUBROUTINE  ITERATE 
COMMON/INPUT/THETA. IN  OMITI98) 

COMMON/WG  INPUT/TYPE  ITER. MAX  DELTA, WGIN  SKIP( 88 ) .THETA  INC* 

$  WGIN  OMIT(2)»DTHETA,WGIN  NO( 5 ) 

COMMON/R/R  SKIP! 18) ,R11,R22*R  0MITI4) 

COMMON/F/F .DFDTHETA.R  BAR11.R  BAR22 

COMPLEX  THETA, FO,F, DFDTHETA, DELTHETA, I ,R11,R22*R  BARI  1 *R  BAR22* 
$  EIGEN.DTHETA 

REAL  MAX  DELTA, lub  REAL.LUB  IMAG 

integer  type  iter.rc  print, r  poly*first, fixed  inc 

DATA ( MAX  I TER*2C  )  » ( LUB  REAL=0 ,0 5 ) * ( LUB  IMAG=0,005) 

DATA( I*(0, 0*1,0) ) 


TYPE  ,13) 
=  ,F5.2) 


11  THETA  =  THETA-DTHETA 
CALL  COMPUTE  F 
FO  =  F 

THETA  =  THETA-fDTHETA 
CALL  COMPUTE  F 
DFDTHETA  =  ( F-FO ) /DTHETA 
DELTHETA  *  -F /DFDTHETA 
FIXED  INC  =  1 
GO  TO  20 

15  FO  s  F 

CALL  COMPUTE  F 
DFDTHETA  «  ( F-FO ) /DELTHETA 
DELTHETA  *  -F /DFDTHETA 
FIXED  INC  *  0 

20  IF ( FIRST  ,EQ,  0)  GO  TO  21 

CALL  PRINT  RS  S  CALL  PRT  THETA  S  FIRST  ■  0 

21  DEL  REAL  *  DELTHETA  $  ABS  REAL  «  ABSF < DEL  REAL) 

IF ( ABS  REAL  ,GT.  THETA  INC)  DELTHETA  *  DELTHETA*THET A  INC/ABS  REAL 
DEL  IMAG  =  -I*DELTHETA  $  ABS  IMAG  =  ABSF(DEL  IMAG) 

I F ( ABS  IMAG  ,GT.  THETA  INC)  DELTHETA  =  DELTHET A*THET A  INC/ABS  IMAG 

THETA  «  THETA+DELTHETA 

CALL  PRT  THETA 

NR  ITER  =  NR  ITER+1 

IF ( NR  ITER  ,GT.  MAX  ITER)  GO  TO  30 

DEL  REAL  =  DELTHETA  $  DEL  IMAG  ■  -I*DELTHETA 

ABS  REAL  =«  ABSF  (DEL  REAL)  $  ABS  IMAG  =  ABSF  ( DEL  IMAG) 

IF ( ABS  REAL  ,LE,  LUB  REAL  ,AND.  ABS  IMAG  ,LE •  LUB  IMAG)  GO  TO  30 
IF ( ABS  REAL  ,LE •  MAX  DELTA  .AND.  ABS  IMAG  ,LE.  MAX  DELTA)  15,11 

30  IF ( FIXED  INC  ,E0.  1)  GO  TO  31 
THETA  =  THETA-DTHETA 
CALL  COMPUTE  F 
FO  =  F 

THETA  =  THETA-fDTHETA 


PRINT  101, TYPE  ITER 

101  F0RMAT(1H1,14HITERATI0N 
PRINT  102.MAX  DELTA 

102  FORMAT ( 1H  , 1 2HMAX  DELTA 
FIRST  =  1 

NR  ITER  =  0 
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CALL  COMPUTE  F 
DFOTHETA  =  ( F-FO ) /DTHETA 

31  IF ( TYPE  ITER  *EQ.  1)  DFOTHETA  «  CR  BAR22*R22-1*0 )*DFDTHETA 
IF ( TYPE  ITER  .EQ.  2)  DFDTHETA  «  tR  BARI 1«R1 1-1*0 )*DFDTHETA 


501 

502 

503 


CALL  PRINT  RS 
PRINT  501 #NR  ITER 
FORMAT ( 1H0 .23HI TERAT IONS 
PRINT  502.DELTHETA 
FORMAT ( 1H0  » 1 1HDEL THETA  ■ 
PRINT  503 #F 

FORMAT ( 1H0 • 1 3HF  FUNCTION 
RETURN 


PERFORMED  «  *I2» 

*C ( FI 1 *7*F1 1*7 1) 

*  *C(F13*5*F13*5) ) 


END 
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SUBROUTINE  COMPUTE  F 

COMMON /WG  INPUT/TYPE  ITER tWGI  N  SKIP<6)*R  POLY.WGIN  OMITI92) 

COMMON/R/R  SKIP( 18) »R11*R22*R!2*R21 

COMMON /F /F  #F  SKIP(2)*R  BAR11.R  BAR22 

COMPLEX  R11.R22.R12.R21.F.R  BAR11.R  BAR22 

INTEGER  TYPE  I  TER  *R  POLY 


CALL  R  BARS 

IF ( R  POLY  .EO.  0)  CALL  INTEG 
IF ( R  POLY  .EQ.  1)  CALL  USE  POLY 

IFITYPE  ITER  .EO.  0)  F  ■  (R  BAR11#R11-1.0)*(R  BAR22*R22-1 .0 ) 
$  -R  BAR11*R  BAR22*R12*R21 

IF ( TYPE  ITER  .EO.  1)  F  «  R  BAR11*R11-1 .0 
IF ( TYPE  ITER  .EQ.  2)  F  •  R  BAR22*R22-1 .0 
RETURN 

END 
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subroutine  r  bars 


COMMON/INPUT/THETA, FREQ  KHZ.IN  SK IP < 3 ) * ALPHA*H,D» IN  0MITI91) 

COMMON /WG  INPUT/WG  SKI P ( 2 ) ,EPS I  LON , EPSI LON  0»SIGMA,WG  0MITI95) 
COMMON/F/F  SK I P { 4 ) *R  BARI 1 ,R  BAR22 
COMMON/RB  CP/S  H  FACTOR 

COMPLEX  THETA, S  H  FACTOR* I »NG  SQ*C,S,SQ  ROOT .RATIO  RT, 

$  PD *QD *H1  D,H2  D*H1  PRIME  D,H2  PRIME  D,CAP  HI  D,CAP  H2  D* 

$  PO *QO , HI  O.H2  0.H1  PRIME  0,H2  PRIME  0,CAP  HI  OtCAP  H2  0, 

S  A1ST  *A2ND*A1*A?»A3*A4*R  BAR  1 1 *R  BAR22  * 

S  C  POWER, Z1»Z2,C  TEMP, 

$  CPLX  sqrt.cplx  expf.cplx  SINF.CPLX  COSF 

REAL  K.LOG  K  OVR  A*K  OVER  A  OT,K  OVER  a  TT,ND  SQtNO  SQ 
DATA ( PI *3,141592653) 

DATA{ VEL  LIGHT=2.997928E05) 

DATA! I*<0.0,1.0 I ) 

DATA ( DEG  TO  RAD*0. 01745329252 ) 

DATAITEST  TH  IM=10.0) 


OMEGA  *  2 • 0*P I*FREQ  KHZ*1000.0 
K  «  OMEGA/VEL  LIGHT 

NG  SQ  »  ( EPSI  LON— I*SIGMA/OMEGA ) /EPS I  LON  0 

C  «  CPLX  COSF ( THETA*DEG  TO  RAD) 

S  *  CPLX  SINF ( THETA*DEG  TO  RAD) 

LOG  K  OVR  A  «  LOGF(K/ ALPHA) 

K  OVER  A  OT  *  EXPF ( LOG  K  OVR  A/3,0) 

K  OVER  A  TT  »  K  OVER  A  OT»*2 

A  OVER  K  OT  *  l.O/K  OVER  A  OT 

A  OVER  K  TT  *  A  OVER  K  OT**2 

ND  SQ  «  1. 0-ALPHA* (H-D) 

NO  SQ  «  1 ,0-ALPHA*H 

SQ  ROOT  «  -CPLX  SORTING  SQ-S**2) 

RATIO  RT  «  NO  SQ/NG  SQ*SQ  ROOT 

PD  «  K  OVER  A  TT*(1,0-ALPHA*(H-D)-S**2» 

PO  ■  K  OVER  A  TT*( 1,0-ALPHA*H-S**2) 

0  TERM  ■  K  OVER  A  TT* ( ALPHA/K ) *»2 /4 . 0 
OD  *  PD-Q  TERM 
00  »  PO-Q  TERM 

THETA  IM  ■  -I*THETA 

IF (-THETA  IM  ,GT,  TEST  TH  IM  .OR,  D  #EQ,  0,0)  GO  TO  50 
CALL  MD  HANKEL ( QD ,H1  D,H2  D,H1  PRIME  D,H2  PRIME  D) 

CALL  MD  HANKEL 1 00, HI  0»H2  0,H1  PRIME  0*H2  PRIME  0) 

CAP  HI  D  »  HI  PRIME  D+A  OVER  K  TT*H1  D/2.0 

CAP  HI  0  ■  HI  PRIME  O-fA  OVER  K  TT*H1  0/2,0 

CAP  H2  D  «  H2  PRIME  D+A  OVER  K  TT*H2  D/2.0 

CAP  H2  0  ■  H2  PRIME  0+A  OVER  K  TT*H2  0/2.0 

A1ST  «  CAP  H2  0-I*RAT 10  RT»K  OVER  A  0T*H2  0 
A2ND  *  CAP  HI  0-I*RAT 10  RT»K  OVER  A  0T*H1  0 
A1  «  C*ND  SQ*(-H1  D*A1ST+H2  D*A2ND) 

S  H  FACTOR  *  I-A1ST*H1  0+A2ND*H2  0)/(-AlST*Hl  D+A2ND*H2  D) 


1*0 


1 


A1ST  =  - I*A  OVER  K  OT*CAP  H2  O-RATIO  RT*H2  0 
A2NO  =  -I*A  OVER  K  OT*CAP  HI  O-RATIO  RT*H1  0 
A2  =  -CAP  HI  D*A1ST+CAP  H2  D*A2ND 
R  BARI  1  *  (A1-A2)/(A1+A2) 


CALL  MD  HANKELIPD.Hl  D*H2  D*H1  PRIME  D.H2  PRIME  D) 

CALL  MD  HANKELIPO.Hl  0*H2  0*H1  PRIME  0»H2  PRIME  0) 

A1ST  =  I *A  OVER  1C  0T*H2  PRIME  O+SQ  R00T*H2  0 

A2ND  =  I*A  OVER  K  0T*H1  PRIME  O+SQ  R00T*H1  0 

A3  =  -HI  PRIME  D+A1ST+H2  PRIME  D+A2ND 
A1ST  «  H2  PRIME  0-I*K  OVER  A  OT*SQ  R00T*H2  0 

A2ND  *  HI  PRIME  0-I*K  OVER  A  OT*SQ  R00T*H1  0 

A4  *  C*(-H1  D*A1ST+H2  D*A2ND ) 

R  BAR22  =  ( A3+A4 ) / ( A4-A3 ) 

RETURN 

C  FLAT  EARTH 

50  C  POWER  «  CPLX  EXPF ( -2 .0* I *K*D*C ) 

R  BARI  1  «  (NG  SQ*C-SQ  R00T)/(NG  SQ*C+SQ  ROOT)*C  POWER 
R  BAR22  *  (C-SO  ROOT)/(C+SQ  ROOT)*C  POWER 
Z1  *  C+SQ  ROOT/NG  SQ 

22  *  -C+SQ  ROOT/NG  SO  , 

C  POWER  «  CPLX  EXPF ( I *IC*D*C ) 

C  TEMP  «  ( Z1-Z2 ) / ( Z1*C  P0WER-Z2/C  POWER  I 

S  H  FACTOR  =  C  TEMP**2 

RETURN 


END 


J»1 


subroutine  comp  proc 


COMMON/ INPUT /THETA* FREQ  KHZ. IN  SK IP( 3 )  » ALPHA *H *D , I N  0MIT(91) 
COMMON/WG  INPUT/WGIN  SK I P ( 3 ) .EPS  I  LON  O.WGIN  OMIT(2>* 

S  REFL  HT.WGIN  NO(93) 

COMMON/F/F  SKIP<2) .DFDTHETA.R  BAR  1 1 *R  BAR22 
COMMON/R/R  SKIP! 18) *R11 .R22.R 12 »R21 
COMMON /RB  CP/S  H  FACTOR 

COMPLEX  I. S. SORT  S.E  I  PI  OV  4. STORE  1. STORE  2. 

S  S  H  FACTOR. SCRIPT  H, 

$  DFDTHETA. THETA. CPLX  SQRT.CPLX  SINF.S  THETA  P.C  THETA  P. 

S  R11.R22.R12.R21 »R  BARI  1  .R  BAR22 • 

S  EXTRA, WAITS  EF ,DENOM .POLAR 

REAL  K.LOGEIO, I MAG 
DATA! I*(0. 0,1.0) ) 

DATA! DEG  TO  RAD=0. 01745329252 )•( RAD  TO  DEG=57. 2957795 1 ) 

DATA! VEL  LIGHT*2.997928E05) 

DATA(PI«3. 141592653) 

DATA! E  I  PI  OV  4* ( 0. 707106781 2 ,0 • 7071 067812 ) ) 

DATA (L0GE10*2. 302585093) 


K  «  2.0*PI*FREQ  KHZ* 1000.0 /VEL  LIGHT 
CAP  K  *  1 • 0/  ( 1 • 0-ALPHA#H/2.0 ) 

S  *  CPLX  SINF(THETA#DEG  TO  RAD) 

S  THETA  P  *  CAP  K*S 

STP  REAL  «  S  THETA  P 

STP  IMAG  «  -I*S  THETA  P 

PHASE  VEL  «  VEL  LIGHT/STP  REAL 

V  OVER  C  *  PHASE  VEL/VEL  LIGHT 

ATTEN  *  -8.6858896*K#STP  IMAG*1000.0 

SCRIPT  H  ■  EXPF(-ALPHA*D/2.0)*S  H  FACTOR 
SORT  S  «  -CPLX  SORT  (  S ) 

STORE  1  «  K*K*SORTF(K)*S*SQRT  S*E  I  PI  OV  4 
$  / ( 2.0*EPSI LON  0*SORTF(2.0*PI )  ) 

STORE  2  *  (1.0+R  BAR  11 )*#2*(  1.0-R  BAR22*R22 ) *SCR I PT  H**2 
$  / { DFDTHETA /DEG  TO  RAD*R  BAR11 ) 

EXTRA  «  STORE  2*S 

WAITS  EF  «  -I *K*REFL  HT/2.0*EXTRA 

CALL  MAG  ANGLE! EXTRA, EXTRA  MAG, EXTRA  ANG) 

EXTRA  ANG  *  EXTRA  ANG*DEG  TO  RAD 

CALL  MAG  ANGLE (WAITS  EF, WAITS  MAG, WAITS  ANG) 

WAITS  ANG  «  WAITS  ANG*DEG  TO  RAD 

WAIT  EF  DB  »  20.0*LOGF (WAITS  MAG ) /LOGEIO 

DENOM  »  R  BAR11*R11-1.0 
REAL  »  DENOM 
IMAG  «  -I*DENOM 

IF( REAL  .EO.  0.0  .AND.  IMAG  .EO.  0.0)  DENOM  ■  1.0E-150 

POLAR  «  -R  BAR 11*R 12 /DENOM 

CALL  MAG  ANGLE! POLAR ,POLAR  MAG, POLAR  ANG) 

C  THETA  P  «  CPLX  SQRT(1.0-S  THETA  P**2) 

CTP  REAL  ■  C  THETA  P 

COS  «  SOF<TF  (CTP  REAL**2-STP  IMAG**2) 

THETA  P  RL  ■  ACOSF ( COS )*RAD  TO  DEG 


k2 


i 
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THETA  P  IM  =  LOGFUCTP  REAL+STP  I  MAG ) /COS ) *RAD  TO  DEG 


PRINT  100 

100  FORMAT ( 1H0,//// ) 

PRINT  101 *PHASE  VEL 

101  FORMAT < 1H  ,2X,17HPHASE  VELOCITY  »  ,E12.5,11H  KM  PER  SEC) 

PRINT  102 *V  OVER  C 

102  FORMAT ( 1H  ,2X,24HPHASE  VELOCITY  OVER  C  =  ,F9.5) 

PRINT  103, ATTcN 

103  FORMAT ( 1H0  »2X , 14HATTENUAT ION  =  .E12.5*3H  DB) 

PRINT  104, EXTRA  MAG 

104  FORMAT ( 1H0,2X» 12HEXTRA  MAG  =  .E12.5) 

PRINT  105*EXTRA  ANG 

105  FORMAT ( 1H  »2X«14HEXTRA  ANGLE  =  .F10.5.4H  RAD) 

PRINT  112* WAITS  ANG 

112  FORMAT ( 1H0 *2X * 35HPHASE  OF  WAITS  EXCITATION  FACTOR  *  tF!0.5#4H  RAD) 
PRINT  113 «WA I T  EF  DB 

113  FORMAT ( 1H  .2Xt32HWAITS  EXCITATION  FACTOR  IN  DB  =  *E12.5t3H  DB ) 
PRINT  116.P0LAR  MAG 

116  FORMAT ( 1H0  #2X  * 19HP0LAR I Z AT  ION  MAG  =  *E12.5) 

PRINT  117, POLAR  ANG 

117  FORMAT ( 1H  ,2X,21HPOLARIZATI0N  ANGLE  =  ,F10.5,4H  DEG) 

PRINT  118, THETA  P  RL, THETA  P  IM 

118  FORMAT (1H0,2X,14HTHETA  PRIME  *  ,2F8.3,4H  DEG) 

PRINT  119 

119  FORMAT ( 1H1 ) 

RETURN 

END 
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SUBROUTINE  WG  OUTPUT 

COMMON/INPUT /THETA.THETA  IMG* IN  OMIT (98) 

COMMON/R/R  SK I P ( 18 ) *R ( 4 ) 

COMMON/F/F.DFDTHETA.R  BAR (2) 

COMPLEX  R  BAR.R.ALL  RS. I *F .DFDTHETA 
REAL  IMAG  PART. MAG 
DATA( I=(0.0*1.0) ) 

DIMENSION  ALL  RS( 8) »MAG( 8 ) .ANGLE ( 8 ) .REAL  PART(8).IMAG  PART<8) 


ENTRY  PRINT 
DO  11  J=1 .2 

11  ALL  RS ( J )  = 
DO  12  J-1.4 

12  ALl  RS( J+2 ) 
ALL  RS<?)  * 
ALL  RS<8)  - 


RS 

R  RAR(J) 

=  R  ( J) 

R  BAR ( 1 ) *R { 1 ) 
R  BAR { 2) *R ( 2 ) 


DO  13  J*1 .8 

CALL  MAG  ANGLE ( ALL  RS( J ) .MAG( J ) .ANGLE { J ) ) 
REAL  P ART ( J )  *  ALL  RS(J) 

13  IMAG  PART ( J )  «  -I*ALL  RS(J) 


PRINT  200 

200  FORMAT (1H0.8X.5HTHETA.5X.9H11R11  BAR.6X.7H1R1  BAR .8X.5H11R11. 

$  10X.3H1R1.9X.4H1R11.9X  *4HI 1R 1 .2X . 1 1HR 1 1  BAR*R 11. 

$  4X.9HR1  BAR*R1 / ) 

PRINT  201. THETA. (REAL  PART ( J ) . J*1 .8  ) 

201  FORMAT ( 1H  .4HREAL.F9.3.1X.8E13.5 ) 

PRINT  202. THETA  IMG. (IMAG  PART { J ) . J* 1 .8 ) 

202  FORMAT ( 1H  .4HI MAG.F9. 3 . IX .8E1 3 .5 ) 

PRINT  203 • ( MAG (J)  . J- 1 .8 ) 

203  FORMAT ( 1H  .3HMAG .  1 IX .8F13 • 5 ) 

PRINT  204 . ( ANGLE ( J ) « J*1 .8 ) 

204  FORMAT ( 1H  .5HANGLE.  9X.8F13.5) 

RETURN 

ENTRY  PRT  THETA 

CALL  MAG  ANGLE (F.F  MAG.F  ANGLE) 

CALL  MAG  ANGLE!  DFDTHETA. DFDT  MAG.DFDT  ANGL) 

PRINT  500. THETA. THETA  IMG.F  MAG.DFDT  MAG 
500  FORMAT ( 1H0.8HTHETA  *  .2F 10. 3 . 10X . 8HF  MAG  =  .E10.3.5X. 

$  15HDFDTHETA  MAG  «  .E10.3) 

RETURN 

END 


SUBROUTINE  MO  HANKEL < Z .HI ,H2  .HI  PRIME. H2  PRIME) 


COMPLEX  Z.I.H1.H2.H1  PR1ME.H2  PR  I  ME .ZCUBED  .Z  POWER ♦ TERM1 * TERM2  • 

S  TERM3 .TERM4 .TERM5 .TERM6 .SUM  1 .SUM 2  *SUM3  »SUM4 . SUMS .SUM6 » 

$  SUM7.SUM8.F.G.F  PRIME.G  PRIME. ONE  RT  Z.CPLX  SORT. 

*  Z  FOURTH.SQRT  Z  CUB.  EXP1 . EXP2 .EXP3 . EXP4. EXP5 « 

$  CPLX  EXPF.BETA.RT  Z ,GM2'r  .GPMFP •  I  POWER. M  I  POWER* 

S  Z  M3HALVS.Z  M3HAFf,  M 

DIMENSION  A(23) .8(23) * C ( 2 3 ) »D( 23) .CAP (14) 


DATA! A*  0.9304 

3671 

693. 

31.0145 

5723 

097. 

206.7637 

1487 

316 

574.3436 

5242 

545  . 

870.2176 

5519 

008, 

828.7787 

1922 

864 

541.6854 

3740 

434. 

257.9454 

4638 

302. 

93.4584 

9506 

631 

26.6263 

5187 

0  74. 

6.1210 

0043 

0056. 

1.1592 

8038 

4480 

0.1840 

1275 

9441  . 

0.0248 

3303 

0964, 

0.0028 

8420 

8010 

0.0002 

9133 

4142. 

0.0000 

2582 

7495. 

0.0000 

02C2 

5686 

0.0000 

0014 

1557. 

0.0000 

0000 

8870. 

0.0000 

0000 

0501 

0.0000 

0000 

0026. 

0.0000 

0000 

0001  ) 

DATA(B*  0.6782 

9872 

514, 

11.3049 

7875 

240, 

53.8332 

3215 

431 

119.6294 

0478 

735. 

153.3710 

3177 

865. 

127.8091 

9314 

888 

74.7422 

1821 

572, 

32.3559 

3862 

152. 

10.7853 

1287 

384 

2.8532 

5737 

403. 

0.6136 

0373 

6351, 

0.1093 

7678 

98 

0.0164 

2293 

9955, 

0.0021 

0550 

5122, 

0.0002 

3316 

7788 

0.0000 

2252 

8289. 

0.0000 

0191 

5671  , 

0.0000 

0014 

4470 

0.0000 

0000 

9729. 

0.0000 

0000 

0589. 

0.0000 

OCOO 

0032 

0.0000 

0000 

0002. 

o.ooco 

0000 

0000) 

DATA(C*  0.4652 

1835 

846. 

6.2029 

1144 

619, 

25.8454 

6435 

915 

52.2130 

5931 

140, 

62.1584 

0394 

215. 

48.7516 

8936 

639 

27.0842 

7187 

022. 

11.2150 

1940 

796, 

3.5945 

5750 

255 

0.9181 

5006 

451  , 

0.1912 

8126 

3439. 

0.0331 

2229 

6699 

0.0048 

4244 

1038. 

0.0006 

0568 

3682, 

O.OOCO 

6555 

C 1 82 

0.0000 

0619 

8599, 

0.0000 

0051 

6550, 

0.0000 

0003 

e220 

0.0000 

COOO 

2528, 

0.0000 

0000 

0150, 

0.0000 

0000 

C008 

0.0000 

0000 

0000, 

0,0000 

0000 

0000) 

DATA(D=  0.6782 

9872 

514, 

45.2199 

1500 

962, 

376.8326 

2508 

015 

1196.2940 

4787 

350. 

1993.8234 

1312 

250. 

2044.9470 

9038 

206 

1420.1021 

4609 

865* 

711,8306 

4967 

351, 

269,6328 

2184 

603 

79.8912 

0647 

290, 

19.C217 

1582 

6880, 

3.7188 

1052 

3339 

0.6076 

4877 

8323, 

0.0842 

2020 

4896. 

0.0100 

2621 

4869 

0.0010 

3630 

1278, 

0.0000 

9386 

7869, 

0.0000 

0751 

2435 

0.0000 

0053 

5074, 

0.0000 

0003 

4135. 

O.OOCO 

OOCO 

1962 

0.0000 

0000 

0102, 

0.0000 

0000 

0005) 

DATA (CAP*0. 1041 

6666 

6666 

6666  7,0,0835 

5034  7222  2222 

2, 

0.1282 

2657 

4556 

3271  6,0. 

2918 

4 9 02  6464  1404 

6  , 

0 . 8d 1 6 

2726 

7443 

7576  5,3. 

3214 

0828  1862  768, 

14.9957  6298  6862 

6,78.9230  1301  1587 

,474.4515 

3886 

e. 

$  3207.4900  91.2  4086.5496*19  8923.12*179  1902. C. 

$  1748  4377.0) 


DAT  A ( PI «3. 1415  92653  58979) 

DATA( I»(0. 0.1.0 ) ) 

DATA ( ALPHA=0.85 3  667  ?18  838  951  ) 


^5 


CALL  MAG  ANGLE! Z *ZMAG*PHI ) 

2  POWER  *  1.0 
Z  CUBED  *  2**3 

IF(ZMAG  .GT.  4.2 )  GO  TO  50 

IF ( ZMAG  .LT.  3.2)  21*22 

21  N*12 

GO  TO  30 

22  IF ( ZMAG  .LT.  4.1)  23*24 

23  N*15 

GO  TO  30 

24  N*23 

30  SUM1  «  0.0 
SUM2  *  0.0 
SUM 3  «  0.0 
SUM4  »  0.0 

DO  31  M«1*N 

TERM1  *  A ( M)*Z  POWER 

TERM2  *  B(M)*Z  POWER 

TERM3  «  C  C  M) *2  POWER 

TERM4  «  DIM) *2  POWER 

SUM1  *  SUM1+TERM1 

SUM2  -  SUM2+TERM2 

SUM 3  »  SUM3+TERM3 

SUM4  ■  SUM4+TERM4 

Z  POWER  «  -Z  CUBED/200. 0*Z  POWER 

31  CONTINUE 

F  «  SUM1 

G  «  Z*SUM2 

GM2F  «  G-2.0*F 

F  PRIME  *  -Z**02*SUM3 

G  PRIME  «  SUM4 

GPMFP  »  G  PRIME-2. 0*F  PRIME 

RT  THIRD  «  SQRTFIl. 0/3.0) 

HI  «  G+I*RT  THIRD*GM2F 

HI  PRIME  »  G  PRIME*I*RT  THIRD*GPMFP 

H2  «  G-I*RT  THIRD*GM2F 

H2  PRIME  *  G  PRIME-I*RT  THIRD»GPMFP 

RETURN 


50  SUM 5  >1.0 
SUM 6  *  1.0 
SUM?  «  0.0 
SUM 8  *  0.0 
Z  REAL  «  Z 
Z  IMAG  *  -I*Z 
I  POWER  «  I 


M  I  POWER  =  -I 
RT  Z  «  CPLX  SORT ( Z ) 

IF (Z  IMAG  .LT.  O.O)  RT  Z  *  -RT  Z 
ONE  RT  Z  =  1.0/RT  Z 
Z  FOURTH  *  CPLX  SORT (RT  Z) 

REAL  Z  4TH  «  Z  FOURTH 

IF(REAL  Z  4TH  .LT.  0.0)  Z  FOURTH  «  -Z  FOURTH 

BETA  «  ALPHA/Z  FOURTH 

SORT  Z  CUB  *  RT  Z**3 

Z  M3HALVS  =  l.O/SORT  Z  CUB 

Z  H3HAFS  M  =  Z  H3HALVS 

EXP1  «  CPLX  EXPF ( I*2.0/3.0*SQRT  Z  CUB) 

EXP2  «  EXP1*CPLX  EXPF(-I*5.0/12.0*PI  ) 

EXP3  «  1.0/EXP1*CPLX  EXPF( 1*5.0/12.0 *P I ) 

EXP4  «  EXP 1*CPLX  EXPF( 1*11. 0/12. 0*PI  ) 

EXP5  ■  1.0/EXP1*CPLX  EXPF(-I*11.0/12.0*PI ) 

00  51  M* 1*14 

TERM5  «  I  POWER*CAP(M)*Z  H3HAF S  M 
TERM6  «  M  I  POWER*CAP(M)*Z  M3HAFS  M 
SUM 5  «  SUM5+TERM5 
SUM 6  *  SUM6+TERM6 
EM  *  M 

SUM7  ■  SUM7+  (-1 . 5*EM*1 .0/Z )  *TERM5 
SUM 8  «  SUM8-M-1.5*EM*1.0/Z)*TERM6 
Z  M3HAFS  M  «  Z  M3HAFS  M*Z  M3HALVS 
I  POWER  «  I  POWER* I 
M  I  POWER  «  M  I  POWER* ( -I ) 

51  CONTINUE 

IF C Z  REAL  .GE.  0.0  .OR.  Z  IMAG  .GE.  0.0)  61.62 

61  HI  «  BETA*EXP2*SUM6 

HI  PRIME  *  BETA*EXP*(SUM6*(-O.25*1.0/Z+I*RT  Z)+SUM8) 

GO  TO  70 

62  HI  «  BETA* ( EXP2*SUM6*EXP5*SUM5 ) 

HI  PRIME  »  BETA*(EXP2*( SUM6* ( -0. 25* 1 .O/Z+I *RT  Z ) +SUM8 ) +EXP5*( SUM5 
S  *(-0.25*1. 0/Z-I*RT  Z )*SUM7 )  ) 

70  IF (Z  REAL  .6E.  0.0  .OR.  Z  IMAG  .LT.  0*0)  71*72 

71  H2  «  BETA*EXP3*SUM5 

H2  PRIME  »  BETA*EXP3*(SUM5*(-0.25*1.0/Z-I*RT  ZJ+SUM7) 

RFTURN 

72  H2  ■  BETA* (EXP3*SUM5+EXP4*SUM6 ) 

H2  PRIME  «  BETA*{ EXP3* ( SUM5* (-0.25*1.0/Z-I*RT  Z ) +SUM7 J+EXP4* ( SUM6 
S  *(-0.25*1.0/Z+I*RT  Z ) +SUM8 ) ) 

RETURN 

END 


SUBROUTINE  R  POLYNOM 


COMMON/ INPUT /THETA, IN  OMIT (98) 

COMMON/WG  INPUT/WGIN  SKIP(8).T  L I  ST < 1 1 ) ,WGI N  OMIT(70) 
COMMON/R/LOG  R(4),R  SIC  IP  ( 10 )  ,R  ( 4 ) 

COMMON/ 1  NIT  R/ADJ  FLAG 

COMPLEX  THETA, LOG  R*R,LOG  R  MTRX,T  LIST.LAGR  INTP.CPLX  EXPF 
INTEGER  DIMENSN* ADJ  FLAG, DEBUG 
DIMENSION  LOG  R  MTRX(10,4) 


ENTRY  GEN  R  POLY 
PRINT  100 

100  FORMAT ( 1H  ,/////, 2 9HR  MATRIX  INTERPOLATION  POINTS) 
M  »  1 

11  THETA  «  T  LIST(M) 

CALL  INTEG 

DO  12  N«l,4 

12  LOG  R  MTRX (M,N)  «  LOG  R(N) 

IF  ( T  LIST(M-t-l)  ,E0.  0,0)  GO  TO  13 
ADJ  FLAG  *  1 
M  *  M+l  S  GO  TO  11 

13  DIMENSN  «  M 
ADJ  FLAG  «  0 

PRINT  300 

300  FORMAT ( 1H0 «6X«6HT  LIST,14X,5H11R11,17X,3H1R1 , 

S  17X,4H1R11 , 19X.4H1 1 R 1 ) 

DO  31  M*l, DIMENSN 

31  PRINT  301  ,T  LIST (M) , ( LOG  R  MTRX ( M,N ) ,N« 1 ,4 ) 

301  FORMAT ( 1H  •  IX, C ( F6.2 ,F6.2 ) .4 ( 2X,C( F10.4,F10.2 ) ) ) 
RETURN 


ENTRY  USE  POLY 
DO  51  N«1  ,4 

LOG  R (N)  «  LAGR  INTPIT  LIST.LOG  R  MTRX ( 1 *N ) .THETA, DIMENSN ) 
51  R(N)  «  CPLX  EXPF ( LOG  R (N ) ) 

RETURN 

END 


hQ 
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FUNCTION  LAGR  INTPIX.Y.XO.N) 


COMPLEX  X.Y.XO.LAGR  INTP.SUM.PROD 
DIMENSION  X(N).Y(N) 

SUM  *  0.0 

DO  12  J*1.N 

PROD  *1.0 

DO  11  1*1. N 

IF  C I  «EO.  J)  GO  TO  11 

PROD  *  PROD* CXO-X<I))/(X(J)-X( 1)1 

11  CONTINUE 

12  SUM  *  SUM+PROD* Y  ( J ) 

LAGR  I NTP  *  SUM 
END 
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SUBROUTINE  INTEG 
BUDDEN 

COMMON / I NPUT /IN  SKIP! 8 ) .D.REL  PREC.TOP  HT.IN  OMIT(89) 
COMMON/FLG  I NPUT /DEBUG.RC  PRINT *FLG  OMIT(23) 
COMMON/R/LOG  R ( 8 ) *DLOGRDH( 8 ) .HT ,DELH »R  OMIT (8) 

COMMON / 1 /L  *HT  L I ST( 256 ) *OMI T1 < 256 *3 ) 

REAL  LOG  R.LOG  N  LIST .LOGR  ©.LOWEST  HT.LOGEIO 
INTEGER  S  FLAG. DEBUG. RC  PRINT 
DIMENSION  LOGR  Om.DLOGRDH  0(8). 

S  DEL  LOGR  0(8) .DEL  LOGR  1(8).DEL  LOGR  2(8) 

DATA(DELHMIN»0. 03125) 

DATA(LOGE 10*2. 302585093) 


FACTOR  «  EXPFC-REL 
EMAX  *  3.0*FACT0R 
EMIN  *  0. 3*FACT0R 
L  *  1 

HT  *  TOP  HT 
DELH  «  DELHMIN 
SAVE  DELH  *  DELH 
CALL  I NIT  COMP 
CALL  S  MATRIX 
CALL  INITIAL  R 
CALL  R  DERIV 
IF (RC  PRINT  .EQ.  1 
100  FORMAT (1H0.15X. 2 1H 
IF(RC  PRINT  .EQ.  1 
IF ( RC  PRINT  .EQ.  1 


PREC*L0GE10) 


•OR.  DEBUG  «GE.  1)  PRINT  100 
R-MATRIX  INTEGRATION  ) 

•OR.  DEBUG  .GE.  1)  CALL  R  COL  HEAD 
•OR.  DEBUG  .GE.  1)  CALL  PRINT  R 


RUNGE  KUTTA 

20  S  FLAG  •  0 

IF(HT  .EO.  HT  L I ST( L+l ) )  L  «  L+l 
IF (HT-DELH  .GE.  HT  LIST(L*1))  GO  TO  21 
DELH  -  HT-HT  LIST(L*1) 

S  FLAG  »  1 

21  IF (HT-DELH  .LT.  D)  DELH  *  HT-D 

DO  31  1*1.8 

LOGR  0(1)  *  LOG  R ( I ) 

31  DLOGRDH  0(1)  *  DLOGRDH(I) 

TRY  AGAIN 

33  DO  34  1*1.8 

DEL  LOGR  0(1)  «  -DLOGRDH  0(I)#DELH 

34  LOG  R ( I )  *  LOGR  0( I )+0.5*DEL  LOGR  0(1) 

HT  -  HT-0.5*DELH 
CALL  S  MATRIX 
CALL  R  DERIV 

DO  35  1-1.8 

DEL  LOGR  1(1)  -  -DLOGRDH ( I )*DELH 

35  LOG  R ( I )  *  LOGR  0( I )*0.5*DEL  LOGR  1(1) 

CALL  R  DERIV 


V 
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DO  36  1*1.8 

DEL  LOGR  2(1)  «  -DLOGRDHt I )*DELH 

36  LOG  R ( I )  *  LOGR  0<I)+DEL  LOGR  2(1) 

HT  =  HT-DELH+0.5#DELH 
CALL  S  MATRIX 
CALL  R  DERI V 

ERROR  *  0*0 
DO  37  1*1.8 

DEL  LOGR  4  *  ( (-DLOGRDH( I )*DELH+DEL  LOGR  0(I))/2.0 
S  +DEL  LOGR  1(1) +DEL  LOGR  2(I))/3.0 

LOG  R ( I )  *  LOGR  0(I)+DEL  LOGR  4 
ERR  *  ABSF (DEL  LOGR  2(I)-DEL  LOGR  4) 

IF (ERROR  .LT.  ERR)  ERROR  *  ERR 

37  CONTINUE 

IF (ERROR  .GT •  EMAX  .AND.  DELH  .GT.  DELHMIN)  38.40 

38  HT  *  HT+DELH 
DELH  «  DELH/2.0 

IF ( DELH  .LT.  DELHMIN)  DELH  «  DELHMIN 
GO  TO  33 

40  CALL  R  DERIV 

IF ( RC  PRINT  .EQ.  1  .OR.  DEBUG  .GE»  1)  CALL  PRINT  R 

IF (ERROR  .LT.  EMIN)  DELH  *  2.0*DELH 

IF ( S  FLAG  *E0.  1)  DELH  -  SAVE  DELH 

SAVE  DELH  *  DELH 

IF(HT  .GT.  D)  GO  TO  20 

IF(RC  PRINT  .EO.  1)  CALL  RATIO  DIF 
RETURN 


END 
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SUBROUTINE  OERIV  EOU 


COMMON/INPUT /THETA, FREG  KHZ,AZ IMUTH,DIP  ANGLE.MAG  FIELD, ALPHA*H, 

S  IN  OMI T ( 92  I 

COMMON/SP  INPUT/NR  SPEC.COEFF  NU(3),EXP  NU( 3 ) * CHARGE ( 3 )  ,M  RATIOI3) 
COMMON /F LG  INPUT/FLG  SKIP(5)*NEG  ELECT, FLG  0MIT(19) 

COMMON/R/LOG  R 11, LOG  R22.LOG  R 12, LOG  R21, 

S  DLNR11DH,DLNR22DH,DLNR12DH,DLNR21DH, 

S  HT,R  SKIP,R11,R22,R12,R21 

COMMON /M/M 1 1 ,M12,M13,M21,M22,M23,M31,M32,M33,C ,S,CSO 
COMMON /1/L,HT  LI  ST ( 256 ), LOG  N  LIST(256.3) 

COMPLEX  I, THETA ,K  OVER  21, 

S  C,S,CSO,ONE  OVER  C,S  OVER  C, 

S  ILY, I MY , INY ,U,USO,D , 

S  Ml 1 ,M12 ,M13 ,M21 ,M22 ,M23,M31 ,M32 ,M33, 

S  T11,T31,T42,T44, 

S  T12  OVER  C,T14  OVER  C*T32  OVER  C,T34  OVER  C,CT4l, 

S  S11A,D11A,S11R,D11B,S12,D12,S21,D21,S22,D22, 

S  B11,B22,B12,B21,R12R21,C12,C21, 

S  DERIVl 1 *DERIV22*DERIV12*DERIV21» 

S  Rll,R22,R12tR21,LOG  Rll.LOG  R2?,LOG  R12.LOG  R21, 

S  DLNR1 1DH,DLNR22DH,DLNR 12DH,DLNR2 1DH, 

S  CPLX  LOGF.CPLX  EXPF,CPLX  COSF,CPLX  SINF 

REAL  MAG  FIELD, M  RATIO, LOG  N,LOG  N  LIST, 

S  LSOYSO,MSOYSO,NSQYSQ ,LMYSQ,LNYSO,MNYSQ 

DIMENSION  X(5),Y(5) * YSQ( 5 )  ,Z  (  5  ) , COLL  FREQ( 5  )  ,ENt 5 ) , 

S  ILY(5),IMY(5) ,INY(5) , LSQYSQ ( 5 ) ,MSOYSQ( 5 ) , NSOYSO ( 5 ) , 

*  LMYSOC  5 ) ,LNYSQ( 5 ) *MNYSQ( 5 ) 

DATA(PI«3. 141592653) 

DATA (DEG  TO  RAD*0. 01745329252 ) 

DATACCOEFF  Y«1 . 758796E1 1 ), (COEFF  X*3,182357E03) 

DATA( VEL  LIGHT«2.997928E05) 

DATA ( I*(0,0«1,0) ) 


ENTRY  INIT  COMP 

OMEGA  -  2,0*PI*FREQ  KHZ*1000,0 
K  OVER  21  »  (OMEGA/VEL  LIGHT )/( 2,0*1) 

SIN  DIP  «  SINF( DIP  ANGLE*DEG  TO  RAD) 

DIR  CCS  L  ■  SIN  DIP*COSF(AZIMUTH*DEG  TO  RAD) 

DIR  COS  M  «  SIN  DIP*SINF( AZIMUTH*DEG  TO  RAD) 

DIR  COS  N  «  -COSF(DIP  ANGLE*DEG  TO  RAT) 

C  «  CPLX  COSF ( THET A*DEG  TO  RAD) 

CSO  ■  C#C 

S  -  CPLX  S I NF ( THET A*DEG  TO  RAD) 

ONE  OVER  C  «  1,0/C 
S  OVER  C  «  S/C 

DO  11  K*1 ,NR  SPEC 

Y(K)  ■  COEFF  Y#CHARGE(j.  )*MAG  F  I  ELD/ (  OMEGA*M  RATIOtK  )  ) 

YSO(K)  »  Y ( 0**2 

ILYCK)  *  I*DIR  COS  L*Y ( K ) 

IMY(K)  ■  I*DIR  COS  M*Y ( K ) 

INVCK)  *  I*DIR  COS  N*Y ( K ) 

LSOYSO(K)  «  DIR  COS  L**2*YS0(K) 

MSOYSO(K)  «  DIR  COS  M**2*YSQ(K) 

NSOYSO ( K )  ■  DIR  COS  N**2*YSQ(K) 

LMYSO(X)  «  DIR  COS  L*DIR  COS  M*YSO(K) 
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) 


LNYSQ(K)  *  DI  i_*DIR  COS  N*YSQ{K) 

11  MNYSQHO  »  DI.  o  M*OIR  COS  N*YSQ(K) 

RETURN 

ENTRY  S  MATRIX 

Mil  *  0.0  $  M12  =  0.0  S  Ml 3  =  0.0 

M21  *  0.0  S  M22  =  0.0  S  M23  *  0.0 

M31  ■  0.0  S  M32  *  0.0  S  M33  =  0.0 

DO  24  K*1*NR  SPEC 

LOG  N  *  LOG  N  L I  ST  (  L+l  *K  )  ♦(  HT-HT  LlST(L-H)) 

S  /( HT  LIST(U-HT  LIST(L+1>> 

S  *(LOG  N  LI  ST  (L  )  -LOG  N  LISTIL+1.IOI 

EN(K)  *  EXPF ( LOG  N) 

X  C  K  >  -  COEFF  X*(1.0E06#EN(K) ) * CHARGE ' K) **2/ ( OMEGA**2*M  RATIO!*)) 
IF ( NEG  ELECT  .EQ.  1)  XU)  =  -XU) 

COLL  FREQU)  *  COEFF  NU( K ) *EXPF ( EXP  NUC K ) #1000. 0*HT } 

Z  U )  *  COLL  FREQU) /OMEGA 

IF ( NEG  ELECT  .EQ.  1)  ZU)  *  -ZOO 

U  »  1 .0—1*2 (K  1 

USQ  ■  U*U 

D  «  -XU)/(U*(USQ-YSQU)  )  ) 

Mil  *  M11  +  (USQ-LSQYSQU)  )*D 
M12  *  M12+(-INYU)*U-LMYSQU)  )*D 
M13  «  M13-MIMYUUU-LNYSQU)  )*D 
M21  «  M21-MINYU)*U-LMYSQU)  )*D 
M22  »  M22+ (USQ-MSQYSQU ) )*D 
M23  «  M23-M-ILYU)*U-MNYSQU)  )*D 
M31  *  M31* ( -IMY C K ) *U-LNYSQU ) ) *0 
M32  «  M32+ ( ILYU )*U-MNYSQ( K )  )  *D 
24  M33  *  M33-MUSQ-NSQYSQU)  )*D 

CURV  TERM  «  ALPHA* ( H-HT I 
Mil  *  Mll-CURV  TERM 
M22  «  M22-CURV  TERM 
M33  «  M33-CURV  TERM 


D  «  1.O/U.0+M33) 

Til  «  — S*M31*D 

T12  OVER  C  «  S  OVER  C*M32*D 

T14  OVER  C  «  (C+M33*ONE  OVER  C)*D 

T31  «  M23*M31*D-M21 

T32  OVER  C  «  C+!M22-M23*M32*D)*ONE  OVER  C 
T34  OVER  C  »  s  OVER  C*M23*D 
CT41  «  ( 1 .0+M1 1-MI 3*M31*D ) *C 
T42  »  M32*M13*D-M1 2 
T 44  *  — S*M13*D 

S11A  *  T11+T44 
DI 1 A  «  T11-T44 
S11B  «  T14  OVER  C+CT41 

D11B  «  T14  OVER  C-CT41 
S12  «  T12  OVER  C+T42 
D12  ■  T12  OVER  C-T42 

521  «  T34  OVER  C+T31 
D21  «  T34  OVER  C-T31 

522  «  C+T32  OVER  C 
D22  *  C-T32  OVER  C 
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RETURN 


ENTRY 

R  DERI V 

Rll  = 

CPLX 

EXPF ( LOG 

Rll) 

R22  * 

CPLX 

EXPFILOG 

R  22) 

R12  = 

CPLX 

EXPF ( LOG 

R12 ) 

R21  = 

CPLX 

EXPFILOG 

R21 ) 

RL  LOG  R12 

=  LOG  R 12 

IFIRL 

LOG 

R12  .LT.  - 

11.0) 

R12 

o 

. 

o 

H 

RL  LOG  R21 

«  LOG  R21 

IFIRL 

LOG 

R21  .LT.  - 

11.0) 

R21 

H 

o 

. 

o 

Bll  «  Rll* ( D1 1 A— D1 IB ) 

B22  *  R22*D22 

B12  *  R 12*D2 1 

B21  «  R21*S12 

R12R21  *  R 12*R21 

C12  «  R12*S2 1 

C21  «  R21*012 


DERIV11  «  B1 1+B12+B21-S1 1B-S1 1B+  I R12R21*D22+C12+C21-D11 A-Dl  IB)  /Rll 
DERIV22  «  B12+B2H-B22-S22-S22*IR12R21*ID11A-D11BH>B12+B2H‘D22)/R22 
DERI VI 2  *  B11+B12+B22+S11A-S11B-S22«MR11#S12+D12)*IR224-1.0)/R12 
DERIV21  *  B114-B21+B22-S11A-S11B-S22+(R11#D21+S21  )*IR22+1 .0)/R21 


DLNR11DH 

DLNR22DH 

DLNR12DH 

DLNR21DH 

RETURN 


K  OVER  2I»DERIV11 
1C  OVER  2I*DERIV22 
»C  OVER  2I*DERIV12 
K  OVER  2I*DERIV21 


END 


5^ 
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SUBROUTINE  INITIAL  R 


COMMON/ INPUT/ IN  SKIP(5)*MAG  FIELD. IN  0MIT(94) 

COMMON /F LG  INPUT/DEBUG. FlG  OMIT (24) 

COMMON/R/LOG  Rll.LOG  R22.LOG  R12.LOG  R21.R  SKIP(IO). 

S  R11.R22.R12.R21 

COMMON/INI T  R/ADJ  FLAG 

COMMON /M/M 11 *M12  »M13  »M21 *M22 .M23.M31 .M32 .M33 .C .S .CSQ 
COMPLEX  Mil .M12 .M13.M21 .M22.M23.M31.M32.M33.C.S.CSQ. 

S  B3.B2.B1.BO.O.C  TEMP. 

S  D1 1. D13. D31. D33. DELTA. P.T.O 1.02 .PI. P2.T1 .T2. 

S  R11.R22.R12.R21  • 

S  LOG  Rll.LOG  R22.L0G  R12.L0G  R21» 

S  FN  SQ.F  ROOT* 

S  CPLX  SORT.CPLX  LOGF.CONJ 

real  mag  field.mag  0 

INTEGER  DEBUG. ADJ  FLAG 
DIMENSION  DIFF(4) .0(4) .P(  2 ) *T(2) 

DIMENSION  PHASE  R(8).PREV  PHAS  { 8 ) 

DATA(PI*3. 141592653) 

DATA (ADJ  FLAG*0 ) 

EQUIVALENCE (0(1) .01) .(0(2) .02) »(P(1).P1).»P(2)*P2). 

S  (T(1).T1).(T(2) *T2) 

EQUIVALENCE (LOG  Rll. PHASE  R) 


IF (MAG  FIELD  .EQ.  0.0)  GO  TO  50 
B3  «  S*(M13«-M31)/(4.0*(1.0+M33)) 

B2  «  (-(CSQ+M33)*( 1.0+M11 >+M13*M3 1- ( 1 .0+M33 ) * ( CSQ+M22 )+M23*M32 ) 

S  /(6.0*  ( 1.0+M33  )  ) 

Bl  «  S*(M12»M23+M21*M32-(CSQ+M?2)#(M13+M31 ) )/(4.0*(1.0*M33) ) 

BO  *( ( 1.0+M11)*(CSQ+M22)*(CSQ+M33)+M12»M23*M31+M13*M21*M32 
S  -M13* ( CSQ+M22 ) *M31- ( 1 .O+Ml 1 )*M23*M32-M12*M2 1*(  CSQ+M33 ) ) 

S  /(1.0+M33) 

CALL  C  JART IC(B3.B2 .Bl .BO .0. DEBUG ) 

DO  12  N*1 .4 

CALL  MAG  ANGLE (Q(N) .MAG  0. ANGLE  Q) 

IF (ANGLE  Q  .LT.  135.0)  ANGLE  0  »  ANGLE  Q+360.0 
12  DIFF(N)  «  ABSF ( ANGLE  0-315.0) 

DO  15  M«2.4 
DO  15  N*M»4 

IF (DIFF(N)  .GT.  DIFF(M-l))  GO  TO  15 

TEMP  ■  DIFF(N)  $  DIFF(N)  *  DIFF(M-l)  $  DIFF(M-l)  «  TEMP 
C  TEMP  ■  0 ( N )  $  Q(N)  =  O(M-l)  S  Q(M-l)  *  C  TEMP 
15  CONTINUE 

DO  21  N-1.2 

Dll  «  1.0+M11-Q(N)»*2 

D13  «  M13+S*Q(N) 

D31  «  M31+S*0(N> 

D33  ■  1.0+M33-S**2 


DELTA  «  D11»D33-D13*D31 

P(N)  *  (-M12*D33+D13»M32) /DELTA 
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T(N>  *  0 ( N ) *P ( N  ) -S* ( -01 1*M32+M 1 2*D3 1 ) /DELTA 
POYNTING  =  CONJ ( PIN ) ) * (Q{ N )*P ( N ) -S*T ( N) )+Q( N ) 

IF ( POYNTING  .LT.  0.0)  PRINT  201 *Q( N ) ♦ POYNTING 
201  FORMAT! 1H0.8HFOR  Q  =  .C( E11.3.E1 1 .3  )  .11HPOYNTING  = 
21  CONTINUE 

OELTA  *  (Tl*C+Pl)#(C+Q2)-( T2*C+P2 ) * ( C+Ql ) 

Rll  =  i (T1*C-P1 >*(C+Q2 )-( T2*C-P2)*I C+Ql ) ) /DELTA 
R22  =  ( (T1*C+P1  )*(C-Q2)-(T2*C+P2)»(C-Q1) )/DELTA 
R 12  =-2.0*C* ( T 1*P2-T2*P1 ) /DELTA 
R21  =-?.0*C* ( 01-02 ) /DELT  A 


22  LOG  Rll 
LOG  R22 
LOG  R12 
LOG  R21 


CPLX  LOGF(Rll) 
CPLX  LOGF ( R22 ) 
CPLX  LOGF ( R1 2 ) 
CPLX  LOGF ( R21 ) 


IF ( ADJ  FLAG  .EO.  0)  GO  TO  36 
DO  35  N«2.8.2 

33  IF (PHASE  R ( N ) — PREV  PHAS(N)  .LE.  PI)  GO  TO  34 
PHASE  R ( N )  =  PHASE  R(N)-2.0*PI  S  GO  TO  33 

34  IF ( PREV  PHAS ( N ) -PHASE  R(N)  .LE.  PI)  GO  TO  35 
PHASE  R(N)  *  PHASE  R ( N ) +2 .0*P I  S  GO  TO  34 

35  CONTINUE 

36  DO  37  N*2 #8.2 

37  PREV  PHAS(N)  *  PHASE  R(N) 

RETURN 


C  RS  BY  FRESNEL 

50  FN  SO  =  1.0+M11 

F  ROOT  «  -CPLX  SQRT( FN  SO-S*S) 

Rll  «  (FN  SO*C— F  ROOT ) / ( FN  SQ*C+F  ROOT) 
R22  *  (C-F  ROOT ) / ( C+F  ROOT) 

R12  *  0.0  S  R21  *  0.0 
GO  TO  22 


END 


.Ell. 3) 


SUBROUTINE  R  OUTPUT 


COMMON/R/LOG  Rll. Rll  ANGLE.LOG  R22.R22  ANGLE* 

S  LOG  R12*R12  ANGLE.LOG  R21*R21  ANGLE* 

$  R  SKIP(8)*HT.R  OM I T  {  9  ) 

REAL  LOG  RlltLOG  R22.LOG  R12.LOG  R21 
DATAIRAD  TO  DEG*57.2957795l ) 


ENTRY  R  COL  HEAD 
PRINT  100 

100  FORMAT ( 1H0 *8X*2HHT. 10X *5H11R1 1, 1 7X *3H1R1 * 16X.4H1R1 1 *16X*4H11R1 ) 
RETURN 


ENTRY  PRINT  R 

Rll  MAG  *  EXPF ( LOG  Rll) 

R22  MAG  «  EXPF ( LOG  R22) 

R12  MAG  ><  EXPF  ( LOG  R12) 

R21  MAG  «  EXPF ( LOG  R21) 

DEG11  *  Rll  ANGLE*RAD  TO  DEG 
DEG22  »  R22  ANGLE*RAD  TO  DEG 
DEG12  *  R 12  ANGLE*RAD  TO  DEG 
DEG21  -  R21  ANGLE*RAO  TO  DEG 

PRINT  200. HT. Rll  MAG.DEGll.R22  MAG. DEG  22.R12  MAG.DEG12. 
$  R21  MAG.DEG21 

200  FORMAT ( 1H  .F 10. 2 .4 ( F 10.4 . F10.2 > > 

RETURN 


ENTRY  RATIO  DIF 

RATIO  =  R 12  MAG/R22  MAG 

DIFF  *  DEG12-DEG22 

PRINT  301. RATIO 

301  FORMAT ( 1H0 . 19H1R1 1  MAG/1R1  MAG  ■  .E12.5) 
PRINT  3C2.DIFF 

302  FORMAT ( 1H  .23H1R11  ANGLE-1R1  ANGLE  *  .F11.5) 
PRINT  303 

303  FORMAT ( 1H1 ) 

RETURN 

END 
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